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ABSTRACT 


This  report  describes  several  potential  applications  of  Kalman  filters  for  advanced  land  fire 
control  systems.  Two  areas  that  are  especially  important  in  the  moving-target  /  moving-platform 
scenario  are  addressed  in  some  detail:  the  tracking  and  trajectory  prediction  of  multiple 
maneuvering  targets  and  the  prediction  of  gun  pointing  angles  in  the  instant  before  firing.  This 
is  particularly  important  in  the  design  of  a  dynamic  muzzle  reference  system.  The  equations  for 
the  filters  are  developed,  simulations  are  described,  and  some  real  data  is  processed  through  the 
muzzle  angle  prediction  filter.  An  architecture  for  a  complete  advanced  land  fire  control  system 
is  proposed. 


r£sum£ 


Ce  rapport  ddcrit  plusieurs  applications  potentielles  de  filtres  Kalman  pour  les  syst&mes  avanc6s 
de  contrdle  d’artillerie  au  sol.  Deux  secteurs  qui  sont  spdcialement  importants  dans  le  s6nario 
plate-forme  mobile/cible  mobile  sont  abordls  dans  certains  details:  la  prediction  de  la  capture  et 
de  la  trajectoire  de  cibles  manoeuvres  multiples  et  la  prediction  de  1’ angle  de  vis6e  du  canon 
au  moment  pr6c6dant  le  feu.  Ceci  est  particulifcrement  important  dans  la  conception  d’un  systbme 
de  reference  dynamique  de  goulot.  Les  equations  des  filtres  sont  developpees,  des  simulations 
sont  decrites,  et  des  donnees  reelles  sont  traitees  par  le  filtre  de  prediction  de  Tangle  du  canon. 
Une  architecture  d’un  syst&me  avance  complet  de  contrdle  d’artillerie  au  sol  y  est  proposes. 
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EXECUTIVE  SUMMARY 


The  Defence  Research  Establishment  Ottawa  (DREO)  was  asked  to  look  at  potential 
applications  of  Kalman  Filtering  technology  to  various  aspects  of  mobile  land  fire  control 
systems,  with  particular  emphasis  on  the  reduction  of  errors  in  the  moving-target,  moving- 
platform  scenario.  Such  errors  are  exacerbated  by  the  high  dynamic  environment,  for  example, 
when  targets  are  tracked  manually  and  when  shells  are  launched  through  a  flexing  gun  barrel. 

Today’s  direct-fire  support  vehicles  (DFSV)  have  large  guns  that  are  inertially  stabilized, 
day  and  night  sighting  systems,  pulsed  laser  range  finders,  and  digital  fire  control  computers  to 
calculate  lead  angles  that  correct  for  target  range  and  crossing  velocity,  weather  conditions,  etc. 
However  increasing  requirements  for  high  accuracy  at  extended  ranges  while  firing  on  the  move 
have  meant  that  a  new  generation  fire  control  system  will  be  required  to  remove  even  more 
sources  of  error. 

Improvements  in  sensor  technology  and  the  decreasing  cost  and  size  of  high-speed 
computers  has  created  opportunities  for  sophisticated  signal  processing  algorithms  aboard 
tomorrow’s  land  vehicles.  Many  tasks  that  crews  find  extremely  difficult  can  often  be  routinely 
done  in  a  computer.  For  example,  accurately  tracking  a  moving  target  when  one’s  own  vehicle 
is  bouncing  wildly  over  rough  terrain  is  nearly  impossible  for  the  gunner  but  is  a  fairly  simple 
exercise  in  target  tracking  for  an  image  processing  computer.  As  well,  present  ballistic 
computers  use  relatively  simple  look-up  tables  because  of  the  limitations  of  on-board  processing 
power.  Modem  computers  can  compute  and  effect  much  more  accurate  gun  laying  angles  that 
take  into  account  more  up-to-the-millisecond  gun  barrel  sensor  information.  However,  at  present, 
few  land  vehicles  exist  with  such  capabilities. 

Kalman  filtering  is  a  signal  processing  technique  that  has  tremendous  potential  for  next 
generation  fire  control.  It  can  address  most  of  the  remaining  error  sources,  provided  that  they 
can  be  mathematically  modelled  and  that  sensitive  sensors  are  installed  to  measure  such  errors. 
Its  predictive  nature  makes  it  ideally  suited  to  estimate  where  the  target  is  going  to  be  when  the 
shell  arrives,  for  example,  or  when  the  gun  barrel  is  going  to  be  pointed  in  the  correct  direction. 

This  report  will  first  summarize  the  fundamentals  of  the  Kalman  filter  algorithm  and  give 
some  simple  examples  of  how  it  is  used  in  common  military  systems.  Next,  a  series  of  problem 
areas  in  mobile  fire  control  systems  that  can  benefit  from  the  application  of  the  filter  will  be 
described.  Two  areas  that  were  studied  in  some  detail,  namely  target  tracking  and  dynamic 
muzzle  referencing  will  be  the  focus  of  subsequent  chapters.  The  equations  used  to  design  the 
filters  will  be  developed;  a  software  simulator  that  implements  the  filter  is  described;  simulations 
of  tracking  performance  are  demonstrated;  and  real  barrel  flex  data  is  processed  through  a 
candidate  pointing  angle  filter  for  a  predictive  muzzle  reference  system.  Broadening  the  scope 
somewhat,  the  next  chapter  outlines  a  potential  architecture  for  an  Advanced  Development  Model 
(ADM)  that  encompasses  these  subsystems  as  well  as  those  related  to  image  processing,  vehicle 
dynamics,  human  interaction,  etc. 
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1.0  INTRODUCTION  TO  ADVANCED  LAND  FIRE  CONTROL 


1.1  BACKGROUND  OF  THE  RESEARCH  PROJECT 


As  part  of  DLAEEM  Task  139  -  "Kalman  Filter  Design  for  Tank  Fire  Control  Systems," 
the  Defence  Research  Establishment  Ottawa  (DREO)  was  asked  to  look  at  potential  applications 
of  Kalman  filtering  technology  in  various  aspects  of  mobile  land  fire  control  systems,  with 
particular  emphasis  on  the  reduction  of  errors  in  the  moving-target,  moving-platform  scenario. 
Such  errors  are  exacerbated  by  the  high  dynamic  environment,  for  example,  when  targets  are 
tracked  manually  and  when  shells  are  launched  through  a  flexing  gun  barrel. 

The  task  description  sheet  specified  the  following  general  aim:  'To  conduct  a  design  and 
feasibility  study  for  an  MBT  (Main  Battle  Tank)  positioning  and  target  tracking  system."  Among 
the  work  items  were  the  following: 

•  develop  dynamical  models  for  tank/target  motion  and  robust  algorithms  for  Kalman 
filter  tracking,  prediction,  error  estimation,  motion  compensation,  sensor  integration; 

•  develop  a  software  simulation  test-bed  to  evaluate  the  algorithms  against  real  or 
simulated  data  in  a  variety  of  motion  scenarios; 

•  analyze  the  simulations  and  identify  the  effects  of  the  algorithms  on  tracking 
performance;  and 

•  recommend  hardware  and  software  configurations  for  possible  advanced  development 
model  construction. 

The  original  tasking  was  very  broad  by  intent  and  it  was  expected  that  research  directions 
would  shift  as  it  became  apparent  which  were  the  primary  error  sources  that  might  be  solved  by 
these  filtering  techniques.  For  example,  the  extent  of  terrain-induced  muzzle  motion  was 
underestimated  until  discussions  were  held  with  several  tank  gun  experts.  Characterization, 
analysis  and  methods  for  compensation  of  these  motions  subsequently  became  one  of  the  largest 
areas  of  study.  On  the  other  hand,  battlefield  navigation,  as  specifically  mentioned  in  the  aim 
of  the  task  description  sheet,  took  on  a  lower  priority  since  it  did  not  affect  fire  control  accuracy 
per  se,  though  it  is  an  area  where  the  application  of  a  Kalman  filter- based,  integrated  land 
navigation  system  would  enhance  the  overall  survivability  of  the  platform  on  the  battlefield. 

In  fact,  the  project  directions  became  so  broad  that  other  establishments  and  several 
contractors  became  involved.  During  the  course  of  this  task,  an  Inter-Establishment  Working 
Arrangement  (IEWA)  was  initiated  to  bring  together  researchers  from  DREO,  DREV  (Defence 
Research  Establishment  Valcartier),  DRES  (Suffield),  and  DCIEM  (the  Defence  and  Civil 
Institute  for  Environmental  Medicine)  that  were  working  in  various  aspects  of  land  vehicles  and 
fire  control.  As  well,  a  major  development  program  (D6374)  was  initiated  to  carry  on  the  work 
and  deliver  a  field-testable  model  of  an  advanced  fire  control  system  that  could  demonstrate  the 
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techniques  of  this  report  and  those  of  the  other  investigators. 

The  work  summarized  in  this  report  was  earned  out  from  May  1988  to  December  1990 
under  project  no.  PCN  0318E. 


1.2  BUILDING  ON  MODERN  DAY  FIRE  CONTROL  SYSTEMS 


Today’s  direct-fire  support  vehicles  (DFSV)  generally  have  good  fire  control  systems, 
especially  when  compared  to  their  forerunners  of  World  War  n.  Gone  are  the  days  when  the 
drivers  of  turretless  tanks  stopped  to  face  each  other  on  the  battlefield  and  fired  volleys  of  hand- 
aimed  shells.  Modem  DFSV’s  have  large  guns  that  are  stabilized  in  azimuth  and  elevation  (to 
allow  a  reasonably  steady  aimpoint  when  the  vehicle  is  moving),  long-range  (out  to  several 
kilometres)  stabilized  day  and  night  vision  sighting  systems,  and  pulsed  laser  range  finders  to 
measure  target  range  to  within  a  few  meters.  Most  have  digital  fire  control  computers  to 
calculate  lead  angles  and  elevation  offsets  that  correct  for  target  range  and  crossing  velocity,  wind 
speed  and  direction,  air  pressure  ammunition  type  and  temperature,  sight  parallax  and  even  the 
static  droop  of  the  gun  barrel  as  it  warms.  However,  increasing  requirements  for  high  accuracy 
at  extended  ranges,  while  firing  on  the  move,  have  meant  that  a  new  generation  fire  control 
system  will  be  required  to  remove  even  more  sources  of  error. 

Given  the  level  of  sophistication  of  the  equipment  in  a  modem  armoured  vehicle,  one  can 
ask  the  fair  question:  is  there  any  more  that  can  be  done?  The  answer  of  course,  is  yes;  several 
countries  are  rigorously  pursuing  mobile  land  fire  control  research  and  several  NATO  and  TTCP 
panels  and  action  groups  have  been  initiated  to  investigate  outstanding  problem  areas  and  identify 
sources  of  errors  that  are  presently  limiting  fire  control  performance.  TTCP  Panel  W,  Action 
Group  10  (WAG-10)  visited  several  US,  UK  and  Canadian  organizations  pursuing  such  issues 
(and  more  as  their  mandate  referred  to  tank  gun  accuracy  in  general).  The  final  report,  (TTCP- 
WAG10  (1989)),  summarized  the  primary  error  sources  and  recommended  the  creation  of  two 
Key  Technical  Areas,  KTA  5-17  "Advanced  Integrated  Tank  Fire  Control,"  and  KTA  2-1 1  "Sabot 
Separation,"  to  further  investigate  them. 

Among  the  findings  of  WAG- 10  was  that  among  the  major  sources  of  error  in  a  typical 
state-of-the-art  main  battle  tank  engaging  a  moving  target,  while  itself  moving  cross  country,  was 
1)  the  inability  of  the  gun  to  remain  pointed  in  the  proper  direction,  and  2)  the  inability  of  the 
gunner  to  maintain  a  sufficiently  accurate  aimpoint  on  a  maneuvering  target.  Though  the 
magnitudes  of  the  errors  are  fairly  small,  on  the  order  of  1  milliradian  or  so  (or  about  1/20  of 
a  degree  or  3  arc  minutes),  in  the  context  of  the  modem  battlefield,  errors  of  even  these 
magnitudes  can  be  disastrous,  especially  if  the  adversary  has  a  comparable  system.  Modem  high- 
intensity  battles  are  fought  at  ranges  from  a  few  hundred  meters  out  to  a  few  kilometers.  Even 
at  a  nominal  battle  range  of  1200  meters,  a  1  mrad  pointing  error  means  a  miss  distance  (from 
the  desired  aimpoint)  of  1.2  meters  (1200  x  tan  0.001).  This  may  mean  the  difference  between 
disabling  the  target  and  missing  it  entirely  (typical  practice  targets  are  2.3  m  wide). 
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Improvements  in  sensor  technology  and  the  decreasing  cost  and  size  of  high-speed 
computers  have  created  opportunities  for  sophisticated  signal  processing  algorithms  aboard 
tomorrow’s  land  vehicles.  Many  tasks  that  crews  find  extremely  difficult  can  often  be  routinely 
done  in  a  computer.  For  example,  accurately  tracking  a  moving  target  when  one’s  own  vehicle 
is  bouncing  wildly  over  rough  terrain  is  nearly  impossible  for  the  gunner  but  is  a  fairly  simple 
exercise  in  target  tracking  for  an  image  processing  computer.  As  well,  present  ballistic 
computers  use  relatively  simple  look-up  tables  because  of  the  limitations  of  on-board  processing 
power.  Modem  computers  can  compute  and  effect  much  more  accurate  gun  laying  angles  that 
take  into  account  more  up-to-the-milUsecond  gun  barrel  sensor  information.  However,  at  present, 
few  land  vehicles  exist  with  such  capabilities.  It  should  be  mentioned,  on  the  other  hand,  that 
there  are  some  tasks  that  are  relatively  simple  for  humans  but  exceedingly  difficult  for  comouters 
-  target  recognition,  for  example.  Research  into  such  algorithms  is  progressing,  however,  and 
may  one  day  be  feasible  in  fire  control  systems. 

Kalman  filtering  is  a  signal  processing  technique  that  has  tremendous  potential  for  next 
generation  fire  control.  It  can  address  most  of  the  remaining  error  sources,  provided  that  they 
can  be  mathematically  modelled  and  that  sensitive  sensors  are  installed  to  measure  such  errors. 
Its  predictive  nature  makes  it  ideally  suited  to  estimate  v.here  the  target  is  going  to  be  when  the 
shell  arrives,  for  example,  or  when  the  gun  barrel  is  going  to  be  pointed  in  the  correct  direction. 
The  filter,  named  after  Rudy  Kalman  who  developed  it  in  the  early  1960’s,  is  a  method  of 
estimating  the  true  state  of  an  imprecisely  known  dynamical  system  from  noisy  measurements. 
Since  this  describes  almost  every  system  that  one  has  to  deal  with  in  practice,  it  is  not  surprising 
that  the  Kalman  filter  has  seen  tremendous  application  in  many  fields;  everything  from  target 
tracking  to  navigation  to  communications  to  economics.  A  very  good  collection  of  both 
theoretical  and  applications  papers  is  available  in  Sorenson  (1985).  The  filter  can  be  relatively 
computationally  intensive  and,  until  recent  generations  of  microprocessors,  its  real-time 
applications  have  often  required  approximate i.s  or  other  ad-hoc  modifications.  Now  however, 
large  size  filters  are  routinely  implemented  in  real-time. 


J.3  THE  SCOPE  OF  THIS  REPORT 


This  report  will  first  summarize  the  fundamentals  of  the  Kalman  filter  algorithm  and  give 
some  simple  examples  of  how  it  is  used  in  common  military  systems.  Next,  a  series  of  problem 
areas  in  mobile  fire  control  systems  that  can  benefit  from  the  application  of  the  filter  will  be 
described.  Two  areas  that  were  studied  in  some  detail,  namely  target  tracking  and  dynamic 
muzzle  referencing  will  be  the  focus  of  the  next  chapters.  The  equations  used  to  design  the 
filters  will  be  developed;  a  software  simulator  that  implements  the  filter  is  described;  simulations 
of  tracking  performance  are  demonstrated;  and  real  barrel  flex  data  is  processed  through  a 
candidate  pointing  angle  filter  for  a  predictive  muzzle  reference  system.  Broadening  the  scope 
somewhat,  the  next  chapter  outlines  a  potential  architecture  for  an  Advanced  Development  Model 
(ADM)  that  encompasses  these  subsystems  as  well  as  those  related  to  image  processing,  vehicle 
dynamics,  human  interaction,  etc. 
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2.0  KALMAN  FILTERING  FUNDAMENTALS 


i.I  CONCEPTS 


A  Kalman  filter  is  a  minimum  variance  estimator  that  can  be  used  to  estimate,  predict 
and/or  smooth  the  true  state  of  an  imprecisely  known  dynamical  system  from  noisy  measure¬ 
ments.  It  is  usually  implemented  in  software  as  a  digital  filtering  algorithm  but  in  the  past, 
steady-state  Kalman  filter  approximations  were  often  implemented  in  analog  hardware. 

It  is  fairly  straightforward  to  describe  the  operation  of  the  filter  in  words.  First,  a 
mathematical  model  is  derived  that  describes  the  state  of  the  system  to  be  estimated.  This  is 
often  in  the  form  of  a  set  of  differential  equations.  The  filter  also  allows  for  a  noise  process  to 
enter  the  state  equations,  whi  a  models  the  random  or  uncertain  nature  of  their  evolution. 
Another  set  of  equations  that  relates  the  sensor  measurements  to  the  system  states  is  derived  and 
knowledge  of  the  type  of  noise  that  the  sensors  generate  is  also  required.  Finally,  the  filter  must 
be  given  a  rough  initial  estimate  of  the  true  state  and  a  guess  of  the  accuracy,  in  the  form  of  an 
initial  error  covariance  matrix,  of  that  initial  state  estimate. 

The  filter  operates  in  a  recursive  mode.  With  its  initial  estimate  of  the  state,  and  the 
differential  equations  describing  the  evolution  of  the  state,  the  filter  algorithm  can  propagate  the 
state  estimate  forward  to  any  point  in  time.  However,  given  the  noise  and  uncertainty  inherent 
in  the  system  model,  the  estimate  will  start  to  diverge  from  the  true  state.  When  sensor  data  is 
available,  the  filter  removes  as  much  noise  from  it  as  possible  by  extracting  the  maximum 
amount  of  information.  It  then  uses  this  information  to  correct  its  estimate  of  the  state.  At  this 
point,  the  error  in  the  state  estimate  will  have  been  reduced  substantially.  The  filter  then 
propagates  its  new  state  estimate  forward  in  time  until  more  sensor  data  is  available. 


2.2  SOME  POSSIBLE  APPLICATIONS  IN  FIRE  CONTROL 

Kalman  filters  appear  in  many  military  systems  in  various  forms.  Often  they  are  only 
simple  approximations  to  the  full  filter  (because  of  the  number  of  computations  that  must  be 
performed  at  each  time  step)  but  they  still  perform  fairly  well.  A  common  example  of  such  an 
approximation  is  the  a-p  or  a-P-y  target  tracker.  Modem  battlefield  computers  now  have  the 
capacity  to  routinely  implement  the  full  set  of  filter  equations  at  a  very  high  rate  enabling  the 
particular  system  to  achieve  better  accuracy  and  reliability. 

There  are  many  systems  aboard  mobile  land  vehicles  that  can  benefit  from  the  integration 
of  a  Kalman  filter  into  the  fire  control  system.  Two  areas  that  have  seen  the  near  universal 
application  of  the  filter  are  multi-sensor  integrated  navigation  systems  and  advanced  target 
tracking  systems.  In  fact,  the  filter  was  developed  for  navigation  systems  in  the  1960’s  and 
nearly  all  integrated  navigation  systems  in  service  today  have  a  Kalman  filter  at  the  heart  of  their 
processing  algorithms.  Tracking  systems  are  ideally  suited  to  the  algorithm  because  of  the 
inherent  filtering  and  prediction  operations  that  are  necessary  for  accurate  tracking  without  loss 
of  lock  during  target  maneuvers  or  obscurations.  The  tracking  application  of  the  filter  will  be 
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expanded  upon  in  a  later  chapter.  Dynamic  muzzle  referencing  is  a  relatively  new  idea  that  is 

also  an  excellent  candidate  for  the  Kalman  filter  algorithm.  These  potential  applications  are 

briefly  described  in  the  following  paragraphs. 

a)  Navigation  -  Accurate  positioning  on  the  battlefield  at  night,  in  adverse  weather  or  on 

featureless  terrain  has  proven  to  be  a  tremendous  advantage  in  modem  warfare.  The 
advent  of  GPS  (the  Global  Positioning  System)  has  allowed  individual  vehicles  to 
obtain  very  accurate  position  and  velocity  information  (at  discrete  times)  very 
inexpensively.  The  integration  (via  a  Kalman  filter)  of  GPS  and  an  inexpensive 
inertial  navigation  system  (INS)  greatly  improves  the  accuracy  and  reliability  of  the 
system  because  the  filter  uses  the  GPS  data  to  help  remove  the  errors  from  the  INS 
and  the  INS  provides  continuous  position,  velocity  and  attitude  information  and  can 
navigate  through  outages  or  obscurations  of  the  GPS  satellites.  A  simplified  INS/GPS 
system  is  used  as  an  example  later  in  this  chapter  to  demonstrate  the  idea  behind 
Kalman  filter  based  sensor  integration. 

b)  Target  Tracking  -  Manual  tracking  of  a  distant  target  through  a  high-power  optical  sight 

is  not  easy.  Even  though  the  mirrors  in  the  sight  are  often  inertially  stabilized, 
providing  a  reasonably  steady  image  of  the  target  and  the  background,  the  sight’s 
eyepiece  and  the  gunner’s  head,  for  that  matter,  are  not  One  can  imagine  trying  to 
look  through  the  eyepiece  and  make  minute  adjustments  with  the  hand  controls  to  keep 
the  target  in  the  center  of  the  sight  while  being  jostled  about  by  the  vehicle  as  it  is 
moving  at  high  speed  over  rough,  uneven  terrain.  Trying  to  keep  one’s  head  steady 
and  to  keep  one’s  hands  from  transmitting  motions  into  the  controls  is  exceedingly 
difficult.  An  image  processing  system  attached  to  the  sight  offers  many  possibilities 
to  improve  performance.  Automatic  tracking  algorithms  can  lock  onto  the  target  and 
keep  it  in  the  center  of  the  field  of  view.  Alleviating  the  gunner  of  this  intensive  task 
can  free  him  for  other  functions.  The  tracking  algorithm  can  be  a  Kalman  filter,  or 
perhaps  a  simplified  variation  of  one,  that  uses  a  mathematical  model  of  the  dynamics 
of  the  target  along  with  azimuth  and  elevation  measurements  of  the  target  from  the 
image  processing  system,  to  determine  its  most  probable  position  at  the  next  time  step. 
This  can  be  fed  back  to  the  image  processing  algorithm  to  aid  it  in  determining  what 
portion  of  the  field  of  view  should  be  concentrated  on  to  locate  the  target  in  the  next 
frame,  especially  if  the  target  has  briefly  disappeared  behind  an  obscuration. 

c)  Target  State  Estimation  /  Trajectory  Prediction  -  Most  current  fire  control  systems  assume 

very  simplistic  target  dynamics  when  computing  lead  angles.  For  example,  the  target 
is  often  assumed  to  be  travelling  at  a  constant  horizontal  velocity  with  no  vertical 
velocity.  This  simplifies  the  lead  angle  computations,  since  horizontal  lead  will  be  a 
function  of  the  angular  rate  of  the  turret  (caused  by  the  gunner  rotating  the  sight  to 
keep  the  target  centered)  and  the  range  to  the  target,  while  the  vertical  lead  is  only  a 
function  of  the  range  (plus  the  usual  shell  flight  compensation  terms,  etc.).  A 
trajectory  prediction  filter,  similar  to,  but  somewhat  more  complicated  than,  a  filter 
used  in  the  tracker  in  an  image  processing  system,  can  be  designed  to  predict  with 
much  more  accuracy,  the  location  of  the  target  at  the  expected  time  of  impact  of  the 
shell.  Such  a  filter  can  estimate  target  acceleration,  as  well  as  velocity,  in  all  three 
directions  -  range,  cross  range,  and  vertical.  This  would  greatly  enhance  the  accuracy 
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of  the  fire  control  system  against  highly  maneuverable  targets  such  as  high-speed 
evasive  ground  targets  and,  especially,  helicopters. 

d)  Stabilization  System  Error  Compensation  -  The  gun  and  turret  stabilization  systems  on 

current  land  vehicles  are  good,  but  they  are  not  perfect.  Under  static  conditions,  a 
modem  electric  drive,  or  even  a  later-generation  hydraulic  drive,  can  position  the  gun 
to  within  a  fraction  of  a  milliradian  of  the  desired  position.  However,  under  dynamic 
conditions  the  gun  control  system  will  be  unable  to  remove  all  of  the  morion 
disturbances.  Current  DFSV’s  often  use  a  simple  "fire  inhibit"  scheme  that  prevents 
firing  until  the  gun  is  sufficiently  aligned  to  the  target  If  these  motions  can  be 
reasonably  well  modelled,  a  Kalman  filter  based  compensation  system  could  be 
designed  to  predict  when  the  stabilization  system  will  have  the  gun  pointing  in  the 
proper  direction. 

e)  Dynamic  Muzzle  Referencing  -  The  terrain-induced  flexing  of  the  gun  barrel  will  result  in 

the  muzzle  of  the  barrel  pointing  in  a  different  direction  than  the  gun  control  system 
assumes  it  is.  These  motions  can  be  reasonably  well  modelled  and  hence  can  be  used 
in  a  Kalman  filter  based  compensation  system  to  predict  when  the  muzzle  will  be 
pointing  in  the  desired  direction.  This  is  the  subject  of  Chapter  6.  The  ignition  of  the 
charge  can  be  delayed  for  a  few  milliseconds  so  that  the  shot  exit  time  coincides  with 
a  barrel  flex  zero-crossing  time. 


2.3  KALMAN  FILTER  EQUATIONS 


The  Kalman  filter  equations  can  be  given  in  several  forms,  since  they  have  been  developed 
for  both  continuous  and  discrete  time,  linear  and  nonlinear,  and  time-varying  or  time-invariant 
systems.  Here  we  will  show  the  case  of  the  linear,  discrete,  time-invariant  form.  The  reader  is 
referred  to  e.g.  Gelb  (1974)  for  other  variations  and  more  complete  derivations. 

Let  xk  denote  a  vector  (generally,  bold  face  will  denote  vectors  and  capitals  will  denote 
matrices,  E(*)  denotes  expected  value  and  :=  denotes  a  definition)  of  dimension  n  representing 
the  true  state  of  a  linear,  discrete-time  system  at  time  index  k: 

xk+ 1  =®kxk  +  wk  V-V 

where  <t>£  is  a  known  nxn  matrix  governing  the  transition  of  the  state  from  time  tk  to  tk+J  and 
wk  is  an  n-dimensional,  zero  mean,  Gaussian  random  vector  with  covariance  Qk  that  enters  the 
system  at  tk  and  perturbs  the  state  in  a  random  manner.  This  type  of  model  is  often  assumed  to 
have  arisen  from  a  continuous  time  system  of  the  form 

xc(t)  =  Fxc(t)  +  wc(t)  (2-2) 
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where  xc(t)  -dxjdt  denotes  the  time  derivative  of  the  continuous  time  signal,  xc(t).  Here,  wc(t) 
is  a  continuous  time,  zero  mean,  uncorrelated  Gaussian  stochastic  process  vector  with  power 
spectral  density  matrix  Qc(t), 


E  [wc(t)wc(x)T)  =  Qc(t)8(t-x) 

where  8 (t)  is  the  Dirac  delta  functional  Equation  (2-2)  can  be  discretized  to  obtain  the 
difference  equation  (2-1)  via: 

<bk  = 

vc(t)di  (2.3) 

It  is  assumed  that  the  sensors  of  the  system  provide  information  about  elements  of  the  state 
vector  at  time  tk  but  that  this  may  be  corrupted  by  noise.  This  is  represented  by 

**  =  Hk*k  +  vk  V-V 


where  zk  is  an  m-dimensional  vector  of  measurements,  Hk  is  a  known  mxn  matrix  relating  these 
measurements  to  the  system  states,  and  vk  is  an  m-dimensional  random  noise  sequence  that 
corrupts  the  sensor  data. 

Certain  simplifying  assumptions  are  usually  made  about  the  noise  sequences,  though  they 
are  not  all  necessary  in  the  most  general  Kalman  filter  formulations  and  can  often  be  relaxed. 
The  measurement  noise  is  assumed  to  be  a  zero-mean,  Gaussian  sequence  with  covariance 
and  is  uncorrelated  with  the  initial  state,  the  state  noise,  and  previous  values  of  itself: 


E(vt)  =  0,  E(vtvj)  =  RkStj 

E(v,wJ)  «  0,  E(v^t0r)  =  0 

where  5kj  is  the  Kronecker  delta  function  defmed  as 


(2-5) 


The  state  noise  is  assumed  to  have  a  similar  form: 


E  (wk)  =  0 


E  =  Qkh  (2-6) 

E(w**o>  =  0 


In  order  to  start  the  Kalman  filter,  an  initial  estimate  of  the  state  and  its  covariance,  P,  must 
be  given.  Often,  the  designer  has  some  knowledge  of  the  statistics  of  the  initial  state  vector  so 
the  initial  state  estimate  is  set  to  the  expected  value  of  the  true  state: 


*0  =  E(*o) 


Pq  -  E{(*0-£o)(xo“*o)r} 


(2-7) 


(In  practice,  though,  the  first  few  measurements  are  often  used  to  estimate  the  initial  states.) 

Given  all  these  assumptions,  die  Kalman  filter  will  produce  the  optimal  linear  estimate  of 
the  true  state  at  time  k,  given  all  available  measurements.  Optimal,  here,  is  defined  as  the 
unbiased  linear  estimator  which  has  the  minimum  variance  or  minimum  mean  square  error.  That 
is,  the  state  estimate  is  chosen  from  the  set  of  all  possible  £  vectors  to  satisfy 

E{x  -  £}  =  0 

T  (2-8) 

and,  min  E{(x-*)r(x-*)} 

£ 

Let  Zk  denote  the  set  of  all  measurements  up  to  and  including  time  k: 

2» k  —  {zo>zi> '-->£*}  (2-9) 


One  can  show  that  the  value  of  £  that  satisfies  (2-8)  is  given  by  the  conditional  expectation  of 
x  given  all  the  measurements  up  to  tk.  This  will  be  denoted  by  £k\k: 

*k\k  :=  E(**|Z*) 

The  covariance  matrix  of  the  estimation  error  at  time  tk  given  measurements  up  to  that  time  will 
be  denoted  by  Pk\k: 


pt\ t  ■' 


Similarly,  the  optimal  estimate  of  x  given  only  the  measurements  up  to  tk,j  is 


:=  ^xk\Zk-0 


and  the  corresponding  error  covariance  matrix  is 


Pk\k-\  := 


The  Kalman  iteration  proceeds  with  an  update ,  followed  by  a  propagate.  A  block  diagram 
of  a  general  Kalman  filter  is  shown  in  Fig.  2-1.  For  the  update,  assume  that  estimates  of  the 
state  vector  and  its  covariance  matrix  are  available  at  tk  that  make  use  of  all  measurements  up 
to  tk.jy  i.e.,  and  P^-i  •  When  a  new  measurement,  z*,  arrives,  the  Kalman  filter  updates 
its  estimate  of  the  state  with  the  following: 


Sensor  Measurement 


Fig.  2-1:  The  General  Kalman  Filter  Algorithm 


A  closer  look  at  this  update  equation  shows  that  the  filter  improves  its  estimate  of  the  state  by 
adding  a  correction  term  to  its  previous  estimate.  This  correction  is  formed  by  taking  the 
residual  (the  difference  between  the  sensor  measurement,  zk,  and  the  estimate  of  what  the 
measurement  should  be,  //***|*_i )  and  then  multiplying  this  by  the  Kalman  gain  matrix,  Kk.  The 
gain  matrix  is  chosen  to  satisfy  (2-8)  and  can  be  shown  to  be  given  by 


(2-11) 


The  filter  also  updates  its  estimate  of  the  error  covariance  matrix  at  this  time  via 


pk\k  =  ~  *|*-i 


(2-12) 


From  this  equation,  it  is  reasonable  to  see  (and,  in  fact,  can  be  shown)  that  the  error  covariance 
matrix  will  become  smaller,  since  a  term  is  being  subtracted  from  its  previous  value.  Thus  the 
Kalman  filter  has  reduced  the  uncertainty  in  its  error  estimate  when  it  has  been  updated  with  a 
new  measurement 

The  purpose  of  the  other  stage  of  the  filter,  the  propagate  stage,  is  to  advance  the  current 
estimates  of  the  state  and  the  error  covariance  ahead  in  time  until  the  next  measurement  becomes 
available  at  '*+/•  This  step  is  relatively  simple  as  it  merely  uses  the  transition  matrix  (or 
dynamics  matrix,  as  in  Fig.  2-1)  that  is  assumed  to  describe  the  system  in  (2-1): 


**♦11*  =  ******  |* 

=  *kpk\k*l +  Qk 


(2-13) 


Equations  (2-10)  to  (2-13),  together  with  the  assumptions  (2-5)  and  (2-6)  and  initial  conditions 
(2-7)  form  the  basic  Kalman  filter  algorithm. 


2.4  EXTENDED  KALMAN  FILTER 


As  was  mentioned  before,  there  are  many  variations  of  the  basic  filter.  A  version  of  the 
Extended  Kalman  Filter  (EKF),  will  be  outlined  in  this  section. 
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Consider  the  case  when  the  measurements  are  nonlinear  functions  of  the  state  vector 


zk  =  hk(*k>  +  vk 


(2-14) 


where  h k{')  can  be  a  general  nonlinear  vector  function  of  x*.  Such  a  system  description  is  very 
common  in  target  tracking  systems.  For  example,  since  the  target  is  often  modelled  in  an  x-y-z 
Cartesian  coordinate  frame  and  the  sensor  is  often  a  range-angle-angle  measuring  radar,  a 
nonlinear  coordinate  transformation  is  required  between  the  state  and  measurement  equations. 

It  can  be  shown  that  a  very  close  approximation  to  the  optimal  Kalman  filter  can  be 
obtained  when  this  non-linear  function,  h  k(-) ,  is  linearized  around  the  current  value  of  the  state 
estimate;  that  is,  the  matrix  is  computed  by  differentiating  hk(x)  with  respect  to  the 

vector  x  and  evaluating  the  result  at  the  current  value  of  : 


^*(**|*-i)  := 


dhk(-x), 

dx 


(2-15) 


where  differentiation  of  one  vector  with  respect  to  another  is  defined  in  the  standard  way  as 

ra/wl  ..  Bf‘ 

hrl#  -  stj 

where  the  subscripts  here  denote  the  corresponding  vector  and  matrix  elements.  This  adds  only 
a  small  amount  of  extra  computation  to  the  filter  and  is  fairly  simple  to  implement  The  Kalman 
filter  update  equations  change  slightly:  Equations  (2-10)  to  (2-12)  become: 


11 


The  propagate  stage,  (2-13),  remains  unchanged  for  the  Extended  Kalman  filter: 


2.5  AN  EXAMPLE:  A  SINGLE  AXIS  INTEGRATED  NAVIGATION  SYSTEM 


A  simplified,  one-dimensional  navigation  example  will  be  described  here  to  show  the 
principles  of  designing  a  Kalman  filter  for  a  multi-sensor  integrated  system  that  provides  a  more 
accurate  and  reliable  solution  than  if  the  sensors  were  used  individually.  Modem  navigation 
systems  often  consist  of  several  independent  sensors,  all  giving  measurements  related  to  position, 
velocity  or  heading  that  the  navigator  must  assemble  into  a  reasonable  estimate  of  his  true 
situation.  The  Kalman  filter  has  been  successfully  applied  to  automate  this  process.  A  large 
system  may  have  one  or  more  inertial  navigation  systems  (INS’s),  a  Global  Positioning  System 
(GPS)  receiver,  other  receivers  for  land-based  radio  wave  systems  such  as  Loran-C  and  Omega, 
a  velocity  sensor  (water  speed  log,  air  speed  data,  doppler  velocity,  or  odometer  for  example), 
heading  sensors  (magnetic  compass,  gyrocompass),  height  sensor  (barometric  altimeter,  radar 
altimeter),  or  a  number  of  others.  (See  McMillan  (1987)  for  details  of  such  a  navigation  system 
as  applied  to  arctic  land  vehicles.)  A  full  Kalman  filter  is  quite  complex,  often  estimating  as 
many  as  50  or  60  state  variables  from  10  to  20  measurements. 

Consider  the  much  simplified  problem  of  integrating  two  sensors  (an  INS  and  a  GPS  set) 
to  estimate  position  p  in  one  direction,  say  latitude  (similar  to  Example  4.2-4  of  Gelb  (1974)). 
In  most  navigation  applications,  the  filter  does  not  estimate  the  position,  per  se,  but  rather  the 
error  in  the  position  as  reported  by  the  INS.  This  is  because  these  errors  can  be  very  accurately 
modelled  in  terms  of  linear  stochastic  processes,  whereas  the  position  itself  is  dominated  by  non- 
random  control  input  (from  the  vehicle  operator)  which  cannot  be  so  accurately  modelled.  Also, 
if  the  aiding  sensor  becomes  unavailable,  the  filter  will  still  be  providing  corrected  INS  position 
outputs,  even  though  the  filter  estimates  of  the  errors  will  slowly  degrade  (a  fail-operational 
system). 

Our  simple  model  of  the  single  axis  INS  will  be  an  accelerometer  with  an  output  ait),  that 
is  integrated  to  give  speed  s(t),  and  then  again  to  give  position  pit): 
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P(0  = 

S(t)  =  J [a(x)dt 


However,  the  accelerometer  has  an  unknown  but  constant  bias  error,  b,  that  is  being  twice 
integrated  so  that  the  indicated  position  has  an  error  that  is  growing  as  t2.  The  Kalman  filter  will 
try  to  estimate  this  error  so  that  it  can  be  removed.  The  INS  error  states  modelled  in  the  filar 
are  INS  position  error  bp,  INS  speed  error  bs,  and  INS  accelerometer  error  6a.  The  relations 
between  the  states  are 

J*8p(t)  =  5  p(t)  =  85(1) 
dt 

—8s(t)  =  8S(t)  =  8 a(t)  (=  b,  the  unknown  bias) 
dt 

1.8a(t)  =  bd(t)  =  0 
dt 

which  can  be  written  in  matrix  form  as 


5p(r) 

0  1  0 

[mo' 

5  m 

- 

0  0  1 

s 

to 

8d(t\ 

"o’ 

o 

P, 

ImoJ 

Note  that  this  is  in  the  form  of  the  continuous  time  system  of  Eq.  (2-2)  with  wc  =  0.  This  system 
can  be  discretized,  at  a  sampling  interval  of  AT  =  tk+rtk,  via  (2-3)  to  get  the  state  transition 
matrix  d>k.  If  we  define  F  as 


F  := 


0  1  0 
0  0  1 
0  0  0 


then 


<Dt  =  eF6J  =  /  +  FAT  + 


F2AT2 

2! 


F3AT3 

3! 


+  ... 


1 

A T 

AT2 

II 

2 

0 

1 

AT 

0 

0 

1 

13 


so  the  discrete  system  corresponding  to  (2-1)  is 


Sp  j 
ds 

\?a\k+l 


1  AT 

0  1 
0  0 


AT2 

~T 

AT 

1 


8p 

8s 

8a 


+  0 


or,  by  making  the  appropriate  association  of  notation, 

**♦1  =  *kxk  +  0 


(2-20) 


The  measurement  for  this  example  will  be  the  difference  between  the  indicated  positions 
of  the  GPS  set  and  the  INS: 


zk  =  pk(GPS)  -  pk(INS)  (2-21) 

The  measurement  from  the  INS  consists  of  the  true  position  plus  the  INS  position  error 

pk(INS)  =  p*  +  8pk 


and  the  measurement  from  the  GPS  will  modelled  as  the  true  position  plus  the  GPS  position 
error,  which  we  call  vk: 


pk(QPS)  =  pk  +  v* 


Thus  the  measurement  equation  (2-21)  becomes 

h  =  vk  -  6P* 
=  [-1  0  0) 


Sp 
8s 

L5aJk 


+  v, 


=  Hxk  +  vk 


which  is  in  the  form  required  in  (2-4)  with  H  :=  [-1  0  0].  The  GPS  positioning  error  will  be 
modelled  as  an  uncorrelated  Gaussian  random  noise  sequence  that  satisfies  (2-5)  with  covariance 

^OPS1 


E  {v*v*}  =  Rk 
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All  that  remains  in  the  filter  design  is  to  specify  the  initial  conditions  on  the  state  vector, 
i0|_j ,  and  the  state  covariance  matrix,  P0|_i.  As  is  usual,  the  initial  state  estimate  is  set  to  zero, 
and  the  initial  covariance  estimate  is  chosen  as  some  physically  reasonable  estimate  of  the  error 
variance: 

roi 


c0|-l 


0|-1 


SPo  0 

0  8  Si 


0 

0 


0  0  8 dg 


The  specification  of  the  Kalman  filter  matrices  for  this  simplified  system  is  complete.  Just 
to  complete  the  design,  some  numerical  values  such  as  those  shown  below  can  be  specified: 


aGPS  =  (10° 

SPq  =  (1000  m)2 

8Sq  =  (10  m/s)2 
8dg  =  (1  mis2)2 


That  is  the  essence  of  the  procedure  involved  in  Kalman  filter  design.  There  are  many 
complications  to  this  in  most  real  systems.  For  example,  the  state  process  noise,  wc,  is  seldom 
assumed  to  be  zero  as  this  can  cause  the  filter  gains  to  go  to  zero  and  the  filter  essentially  stops 
operating.  The  numerical  values  of  the  noise  parameters  are  chosen  after  careful  analysis, 
simulation  and  testing  on  real  data  to  achieve  satisfactory  performance.  The  state  dynamics  and 
measurement  matrices  ( <t>k ,  H)  are  generally  much  more  complex  and  various  precautions  must 
be  taken  to  ensure  numerical  stability  and  to  safeguard  against  invalid  measurements.  The 
sequence  of  operations  can  also  be  more  complex  since  all  aiding  measurements  are  seldom 
available  in  such  a  synchronous  fashion. 


3.0  TRACKING  AND  TRAJECTORY  PREDICTION 


3.1  INTRODUCTION 


This  chapter  summarizes  the  design  of  a  simulation  system  that  was  used  to  evaluate 
various  Kalman  tracking  filter  configurations  that  could  enable  land  platforms  to  track  and  engage 
multiple  maneuvering  targets,  either  passively  or  actively.  Typically  a  direct-fire  support  vehicle 
(DFSV)  acquires  and  must  track  multiple  targets  with  a  passive,  2-dimensional  imaging  sensor 
such  as  a  TV  or  infrared  camera.  When  a  particular  target  is  chosen  to  be  engaged,  an  active 
sensor,  such  as  a  laser  range  finder,  radar,  or  laser  radar,  would  be  used  to  obtain  range  and 
possibly  range-rate  information.  In  most  current  generation  vehicles,  a  laser  range  finder  is  used 
to  obtain  one  range  measurement  just  prior  to  firing  that  is  used  to  calculate  the  elevation  and 
azimuth  lead  angles  using  two  decoupled,  second-order,  constant  angular  rate  filters.  Target 
trajectory  prediction  can  be  improved  substantially  if  a  fully  coupled  3-dimensional  Kalman  filter 
tracking  algorithm  is  used  especially  if  range  is  measured  for  even  a  few  scans. 


3.2  TARGET  MODELLING 


This  section  looks  at  some  of  the  typical  coordinate  frames,  the  transformations  between 
them  and  the  structure  of  some  target  dynamics  models  that  can  be  used. 


3.2.1  Coordinate  Frames  and  Transformations 


A  point  in  three  dimensional  space  can  be  specified  in  a  variety  of  ways  but  usually  in 
either  rectangular  or  spherical  coordinates.  As  was  mentioned  Section  2.4,  the  sensor  of  a 
tracking  system  is  usually  most  naturally  represented  in  a  spherical  frame  but  the  target  knows 
nothing  about  the  sensor  frame  and  is  more  naturally  modelled  in  rectangular  coordinates,  since 
it  is  more  likely  to  travel  in  a  straight  line  than  it  is  to  move  along  the  surface  of  a  sphere. 

Consider  the  coordinate  frames  shown  in  Fig.  3-1.  The  principle  rectangular  frame  is 
defined  by  the  Xq-Yq-Zq  axes  which  are  aligned  with  the  local  geographic  North,  East  and 
Down  directions  with  the  origin  at  the  sensor  location.  The  target  at  point  p  can  be  described 
in  rectangular  coordinates  at  the  point  (PxG>  Pyc » PzJ  or  *n  spherical  coordinates  at  (r,b,e)  where 

b  is  the  bearing  angle  from  XG  (North),  e  is  the  elevation  angle  above  the  horizon,  and  r  is  the 
distance  or  range  from  the  sensor  location  at  the  origin.  The  other  rectangular  frame  shown  as 
Xl-Yl-Zl  is  called  the  line-of-sight-to-target  frame  and  is  obtained  by  rotating  the  geographic 
frame,  first  about  the  Zq  axis  through  an  angle  b  and  then  elevating  it  (rotating  about  the  new 
Yg  axis)  by  the  angle  e.  It  has  the  same  origin  as  the  geographic  frame,  although  it  is  shown 
separated  to  reduce  the  clutter  on  the  diagram.  The  line-of-sight  frame  is  called  such  because 
the  x-axis  of  the  frame  always  points  directly  to  the  target. 


Fig.  3-1:  Target  Coordinate  Frames 

Vectors  in  the  geographic  frame  can  be  expressed  in  the  Iine-of-sight  frame  by 


P*L 

cose  cosh  cose  sinh  -sine 

-sinh  cosh  0 

PxG 

P*c 

_ 

Pxc 

pyG 

P*L. 

sine  cosh  sine  sinh  cose 

1 

j 

— i 

N1 

where  Cgl  is  the  Direction  Cosine  Matrix  (DCM)  that  converts  a  vector  expressed  in  G 
coordinates  to  one  expressed  in  L  coordinates.  Similarly  the  inverse  transformation  is 


pxG 

pyG 

cose  cosh  -sinh  sine  cosh 

cose  sinh  cosh  sine  sinh 

P*L 

pyL 

-  cG 
-•  ll 

P*L 

pyL 

» 

-sJne  0  cose 

i _ 

PZl 
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To  simplify  notation  slightly,  let  p^  denote  an  arbitrary  3-dimensional  vector  that  is  expressed 
in  G  coordinates.  Then  we  simply  write 

P(G)  =  c°  pV 
P®  =  Cq  p< O) 

The  bearing  and  elevation  angles  can  be  expressed  in  terms  of  the  geographic  coordinates  as 

r  =  \  PxG  +  Pyg  +  PzG 

b  =  tan ~x(PxgIPyg)  ^  ^ 

e  =  sin  ~\-pZclr) 


and  similarly  the  reverse  is 


PX(J  ~  r  cose  cos* 

PYg  =  r  cose  si nb  ,  (3-4) 


PzG  ~  ~r 

Since  the  target  is  moving  relative  to  the  sensor,  the  line-of-sight  frame  is  rotating  with  respect 
to  the  geographic  frame.  The  rate  of  rotation  of  the  L  frame  with  respect  to  the  G  frame  will 
be  denoted  as  to uq  and  can  be  written  in  L  frame  coordinates  as 


where  «v  denotes  a  unit  vector  in  the  direction  v,  b  is  the  bearing  slewing  rate  and  £  is  the 
elevation  slew  rate.  The  Coriolus  equation  can  be  used  to  transform  velocities  in  a  rotating  frame 
to  a  non-rotating  frame:  the  rate  of  change  of  a  vector  p  in  the  rotating  frame  as  seen  from  the 
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nonrotating  frame  is  known  to  be  ( e.g .  Goldstein  (1950)  p.  133) 


dp  |  f  dp  } 

P  =  +<»LIGxp 

dt  )G  {  dt  JL 


where  the  notation  ( ,)G  denotes  the  coordinate  frame  from  which  the  vector  is  being  viewed. 
Expressing  all  vectors  of  (3-6)  in  the  same,  say  L,  coordinates  and  defining  velocity  as 


yields 


dp ' 
dt 


Now,  since  the  L  frame  is  aligned  with  its  X  axis  towards  p  which  is  a  distance  r  away,  the  target 
position  and  velocity  in  L  coordinates  are  simply 


=  0 


v p .  5 

dt  L  0 


so  that  by  substituting  these  in  the  Coriolus  equation,  we  get  the  velocity  vector  of  p  as  seen 
from  the  non-rotating  G  frame  (though  it  is  still  coordinatized  in  L  coordinates): 


r  ^  -b  sine  r  ^ 

0  +  e  x  0 

PJl  b  cose  L0. 


=  r  b  cose 


The  appropriate  DCM  can  be  used  to  convert  these  velocities  as  seen  in  the  G  frame  from  L  to 
G  coordinates  via 
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Equation  (3-8)  can  be  inverted  so  that  the  angular  rates  are  obtained  from  the  velocities: 


vYJ(r  cose) 
-v7Jr 


(3-9) 


Similarly  the  Coriolus  equation  can  be  used  to  derive  accelerations  as  seen  in  rotating 
coordinate  frames: 


d_ 

dr 


xp  +  (O^qX 
G 


> G 


^LJG^P 


L 


dp_ 

dt 


xp  +  (D^cXfW^Xp) 
G 
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and  then  choosing  the  L  frame  to  coordinitize  all  vectors  and  defining  acceleration  as 


a  := 


dp 2 


dr 


(3-10) 


Jg 


yields 


a2L 


KD 


d2p 

1 ? 


JG 


£p 

K*  JL 


+  2tog;X 


prJL 


da 


^]? xp(t>  + 


(3-11) 


The  angular  acceleration  of  the  rotating  coordinate  frame  can  be  found  by 

,  .  TO 


(  \ 
d^L/G 

IL) 

d 

-bsine 

£ 

(D  r 

*  J 

dt 

bcose 

bcose-bisiae 

so  substituting  the  individual  L  frame  components  of  this  into  (3-11)  yields 
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■'  .  A 
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k 

b  cose 


0 

Au 

(3-12) 


TO 


f  -  r£2  -  rb2cos2e 
2*'hcose  +  rhcose  -  2rb6sine 
-2f&  -  re  -  rb2 cosesine 


and  finally,  the  DCM  can  be  used  to  convert  the  accelerations  as  viewed  from  the  G  frame  to 
coordinates  in  that  frame: 


KG) 


Yg 

G 


az, 


TO 


Again,  Eq.  (3-12)  can  be  inverted  to  find  the  accelerations  as  seen  by  the  tracking  system: 
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f  \ 
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\eJ 


axL  +  +  ri>2cos2e 

dy  t 

— £_  -  2 Lb  +  2bt\ane 


r  cose 


-flz,  y  •  2 

-  -  2—6  -  6  cose  sine 

r  r 


(3-13) 


3.2.2  Target  Dynamical  Models 

The  mathematical  model  used  to  describe  the  motion  of  the  target  is  of  course  dependent 
on  the  type  of  target.  Tracking  an  evasive  target  through  frequent  obscurations  requires  a  much 
different  model  than  that  used  for  a  cooperative  one.  A  good  selection  of  common  target  models 
is  summarized  in  Bogler  (1990),  Chapter  9.  There  are  roughly  4  main  categories  of  generic 
models  that  are  used.  Following  the  nomenclature  of  Fitzgerald  (1981)  these  are: 

1 .  Random  Walk  Velocity  (RWV) 

2.  Exponentially  Correlated  Velocity  (ECV) 

3.  Random  Walk  Acceleration  (RWA) 

4.  Exponentially  Correlated  Acceleration  (ECV) 

Each  of  these  models  could  be  used  in  either  rectangular  or  spherical  coordinates.  Listed  below 
are  the  full  state  dynamics  equations  that  can  be  used  for  each  of  these  models  in  each  coordinate 
frame  (r  -  rectangular,  s  -  spherical). 

Rl.  Random  Walk  Velocity:  rectangular  (RWV(r)): 


PxG 

Vv 

ag 

PyG 

Vy 

XG 

d 

PzG 

% 

dt 

vx 

% 

vY 

*  G 

% 

V 'y 

% 

0  0  0  1  0 
0  0  0  0  1 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 


r  -i 

0 

Pxc 

o 

PYg 

w 

0 

PZG 

V*c 

T 

W°XG 

% 

Warc 

°Z G 

(3-14) 


where  w  are  the  uncorrelated  white  noise  processes  representing  the  acceleration  components. 
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R2.  Exponentially  Correlated  Velocity:  rec 


PxG 
Pro 
d  PZo 
~divXG 


0  0  0  1  0  0 

0  0  0  0  1  0 

0  C  0  0  0  1 

0  0  0  0  0 

0  0  0  0  -p^  0 

0  0  0  0  o  -p. 


(3-15) 


where  the  P  terms  define  the  correlation  times  of  the  Markov  processes  (as  described  in  Gelb 
(1974))  representing  the  acceleration  components. 


R3.  Random  Walk  Acceleration:  rectangular 


~  mi  tiVATfi  i  >1 


0  0  0  1  0  0  0  0  0  ^ 
000010000  Pyg 
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000000001  vzc 
000000000  axc 

000000000  aYc 

0  0  0  0  0  0  0  0  oj  azG 


(3-16) 


where  j  denotes  "jerk"  or  the  rate  of  change  of  acceleration. 
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S3.  Random  Walk  Acceleration:  spherical  (RWACs)’): 
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S4.  Exponentially  Correlated  Acceleration:  spherical  (ECA(s)): 
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3.3  MEASUREMENT  MODELS 


(3-21) 


Target  measurements  are  almost  always  given  in  spherical  coordinates  as  bearing  and 
elevation  angles  (relative  to  some  known  reference  frame)  and,  if  available,  a  range  from  the 
sensor  to  the  target.  We  will  denote  the  measurement  vector  as  z : 
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These  measurements  will  be  assumed  to  consist  of  the  true  values  of  the  quantities  being 
measured  plus  random  white  noise  components,  v,  consistent  with  the  standard  Kalman  filter 
measurement  model  equation,  (2-4): 


zb 

b 

vb 

2e 
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e 

+ 

ve 

2r 

r 

For  a  state  vector  in  spherical  coordinates,  the  map  between  the  measurement  and  state  is  linear 


1  0  0  0  0... 

0  1  0  0  0... 

0  0  10  0  ... 


(3-22) 


but  for  a  state  vector  in  rectangular  coordinates,  the  relationship  is  nonlinear 
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where  the  nonlinear  functions  are  as  in  (3-3): 


hb(x)  =  tan 
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~Pz, 


sp2xg+p2yg  + Pzc 


hr{x)  = 


\pxG  +  pyG  +  Pzc 


(3-23) 


(3-24) 


For  use  in  an  Extended  Kalman  filter,  this  nonlinear  measurement  vector  will  be  linearized 
about  the  state  vector  estimate  £  =  \pXc  PyG  pzG  —  l7”  according  to  (2-15): 
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which,  after  one  has  taken  the  derivatives,  becomes 
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(3-25) 


*  fac  +  fix  a  +  fizc 


1  :=  fac  +  firc 


(3-26) 


3.4  CHOICE  OF  TARGET  MODEL 


There  has  been  a  tremendous  amount  of  literature  arguing  the  various  merits  and 
shortcomings  of  each  of  these  models.  Countless  variations,  combinations  and  extensions  have 
been  proposed,  studied,  simulated,  tested  and  implemented.  A  number  of  adaptive  schemes  that 
switch  between  models  or  combine  outputs  of  parallel  models  also  exist.  Relatively  recent  survey 
papers  include  Woolfson  (1985),  Chang  and  Tabaczynski  (1984)  and  the  texts  of  Blackman 
(1986),  Farina  and  Studer  (1985)  and  others  are  all  invaluable  sources.  One  early  reference, 
Burke  et  ai  (1976),  used  random  walk  acceleration  models  to  track  a  maneuvering  tank  in  both 
geographic  and  line-of-sight  coordinates.  New  refinements,  such  as  Rouhi  and  Farooq  (1989), 
appear  continuously. 

The  choice  of  target  model  can  best  be  done  when  the  specific  scenario  is  known  and  some 
real  data  is  available  to  evaluate  the  various  models.  A  good  tracking  simulator  will  allow 
different  models  to  be  tested  and  evaluated  so  that  a  suitable  one  (or  more)  is  chosen  for  the 
purpose  at  hand.  The  next  two  chapters  describe  a  simulator  that  was  developed  for  potential  fire 
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control  systems.  It  can  implement  any  of  the  above  models  and  it  has  some  interesting  features 
that  may  be  particularly  useful  for  land  fire  control.  Developing  a  flexible  simulation  package, 
rather  than  extensive  model  evaluation,  was  the  primary  goal  of  the  effort. 

There  are  a  few  rules  of  thumb  that  assist  in  the  choice  of  a  model.  In  general,  the 
inclusion  of  acceleration  states  makes  a  filter  more  able  to  track  sudden  maneuvers  but  at  the 
expense  of  a  noisier  estimate  when  the  target  is  not  accelerating.  On  the  other  hand,  the 
exclusion  of  acceleration  states  will  make  the  filter  estimate  more  accurate  when  the  target  is  not 
maneuvering  but  a  sudden  target  acceleration  may  cause  the  filter  to  loose  the  target.  The 
amount  of  noise,  w,  that  is  assumed  to  enter  the  state  model  can  have  a  similar  effect:  the  larger 
its  covariance,  Qk,  the  noisier  but  more  responsive  the  filter,  the  smaller  the  Qk,  the  quieter  but 
more  sluggish  the  filter.  In  general,  one  often  sees  high  Q  uncorrelated  acceleration  filters  for 
tracking  highly  maneuverable  targets  such  as  fighter  aircraft,  and  low  Q  correlated  velocity  filters 
for  tracking  ships  and  submarines.  Land  combat  vehicles  would  fall  somewhere  in  between.  For 
this  study,  moderate  Q  values  were  used  with  both  velocity  and  acceleration  models. 
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4.0  TRACKING  FILTER  SIMULATIONS 


This  chapter  presents  some  results  from  the  simulation  program  that  was  developed  to  test 
the  various  Kalman  filter  based  tracking  algorithms  presented  in  the  previous  chapter.  The 
software  was  written  in  "C"  and  is  very  generic.  The  initial  Kalman  filter  matrices  and 
parameters  are  read  in  from  files  as  is  the  simulated  sensor  data.  Only  a  few  routines  exploit 
certain  features  or  ordering  of  the  state  vector  elements  and  these  were  grouped  together  to 
remind  the  developer  which  routines  might  have  to  be  modified  if  a  different  state  vector 
formulation  was  used. 


4.1  DATA  GENERATION 

4.1.1  Truth  Data  Generation 

A  pre-existing  program  (TGen)  was  used  to  generate  the  truth  data  for  the  target 
trajectories.  Through  a  plain  language  input  text  file,  the  user  specifies  the  target’s  initial 
position,  velocity  and  heading  and  changes  to  these  values  at  any  subsequent  point  in  time.  The 
program  computes  the  deterministic  position,  velocity,  acceleration,  attitude  and  attitude  rates  of 
the  target  and  saves  these  in  a  data  file. 

4.1.2  Sensor  Data  Generation 

A  simple  program  called  SensData  was  developed  that  would  take  the  data  files 
representing  the  true  target  trajectories  and  produce  simulated  bearing,  elevation  and,  if  requested, 
range  measurements  as  might  be  produced  by  a  generic  imaging  sensor  system.  The  program 
allows  pseudorandom  noise  to  be  added  to  the  true  bearing,  elevation  and  ranges,  permits  the 
turning  on  and  off  of  ranging  information  to  simulate  passive  to  active  mode  switching,  and  can 
simulate  target  obscuration  by  removing  sensor  data  for  arbitrary  intervals.  If  there  is  more  than 
one  target,  it  will  merge  the  sensor  data  associated  with  all  targets  into  one  file  and  sort  by  time 
stamp.  It  can  also  simulate  clutter  by  generating  false  targets  at  any  specified  probability  of  false 
alarm. 

Shown  in  Fig.  4-1  below  (and  in  x-y  form  as  Target  1  in  Fig.  5-1)  is  a  simulated  trajectory 
of  a  land  target  moving  at  a  constant  speed  of  15.0  m/s  (  =  54.0  km/hr)  at  a  range  of  initially 
1500m.  It  recedes  from  the  sensor  for  8  seconds,  makes  a  90  degree  right-hand  turn  (at  r=8s, 
duration  5  seconds),  travels  across  the  sensor  field  of  view  for  another  4  seconds,  then  makes 
another  right  hand  turn  (at  f=17s,  duration  5s)  and  travels  toward  the  sensor  for  the  final  3 
seconds.  The  figure  shows  the  true  range,  bearing  and  elevation  to  the  target  as  well  as  the  noisy 
measurements  of  these  quantities  made  by  the  sensor  at  0.1s  intervals.  The  measurements  are 
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Sensor  data  •  Ground  target  at  1 600  m 


Fig.  4-1:  Simulated  sensor  data  for  a  ground  target  at  1600m,  first  receeding,  then  crossing, 

then  closing  on  sensor  location.  Both  truth  and  noisy  data  are  shown. 

where  the  covariance  values  of  the  measurement  noises  are 

(0.00175  rad)2  0  0 

[  v&(*)  ve(*)  vr<*)  ]7  *  0  (0.00175  rad)2  0 

0  0  (3  m) 
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4.2  TRACKING  FILTER  RESULTS 


This  section  shows  some  representative  results  of  typical  simulations.  The  purpose  is  not 
necessarily  to  show  the  best  tracking  filter  for  the  simulated  trajectory  but  rather  to  show  some 
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of  the  features  that  are  available  in  the  simulator.  Again  the  choice  of  the  "best"  filter  depends 
on  the  application.  In  this  chapter  we  assume  that  the  range  information  is  either  available  for 
the  entire  trajectory,  or  unavailable  for  the  entire  trajectory.  The  next  chapter  describes  a  method 
to  handle  the  case  when  range  data  is  intermittently  available  and  presents  some  simulated  results 
of  such  a  situation. 

The  Kalman  filter  program,  called  "tracker",  was  written  in  C.  It  accepts  a  configuration 
file  with  various  control  parameters,  as  well  as  the  continuous  time  state  dynamics  matrix,  F,  and 
power  spectral  density  matrix,  Qc.  It  also  reads  in  the  assumed  level  of  measurement  noise,  i.e., 
the  R  matrix  and  the  names  of  the  sensor  data  file  and  the  truth  data  file.  It  produces  output  files 
for  the  values  of  the  states,  the  state  errors,  and  residuals,  as  well  as  the  square  roots  of  the 
corresponding  diagonals  of  the  covariance  matrix.  The  program  first  does  a  discretization  of  the 
continuous  time  matrices,  and  thereafter  operates  as  a  discrete  EKF. 

4.2.1  Random  Walk  Velocity  (RWV(r))  Model  -  Range  Available 

In  this  case,  range  information  is  available  and  the  RWV(r)  model,  Eq.  (3-14),  is  used. 
This  is  a  3  dimensional  rectangular  coordinate  frame  model.  The  measurement  equation  is  the 
nonlinear  form  of  (3-23)  so  the  extended  Kalman  filter  formulation  is  used.  The  power  spectral 
density  matrix  of  the  state  process  noise  required  for  the  discretization  equations  (2-3)  is 
estimated  in  Appendix  A  where  it  is  shown  to  be 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

Qlt)  = 

c  0  0  0  3.0  0  0 

0  0  0  0  3.0  0 

0  0  0  0  0  0.0 

Shown  in  Fig.  4-2  below  are  the  errors  in  the  state  estimates  of  each  cartesian  position 
coordinate  as  well  as  the  filter  estimates  of  the  square  root  of  the  covariance  of  that  estimate,  i.e., 

xfk/k)  =  Xj(k)  -  Jtj(kjk)  and  sjP^kJk)  where  the  subscript  refers  to  an  individual  element.  As 
well,  the  velocity  state  estimation  errors  and  covariance  elements  are  shown  in  Fig.  4-3.  The 
measurement  residuals,  Vj(k/k)  =  zfjk)  -  hftik/k))  as  well  as  the  square  root  of  the  filter 

estimated  residual  covariance,  { [HifyPik/tyHik)7]^  +  /?„•(£)} 1/2 ,  are  shown  in  Fig.  4-4. 

The  state  estimate  errors  for  the  most  part  lie  within  their  expected  standard  deviation, 
indicating  a  reasonably  well- tuned  filter.  The  target  right  hand  turns  during  r=8  to  r=13  and  r=17 
to  r= 22  cause  larger  state  errors  to  occur  during  the  maneuvers  and  they  disappear  after  the 
maneuver  is  completed.  The  smaller  Q  value  for  the  vertical  velocity  state  shows  up  in  the 
figures  as  smaller  state  estimation  errors  and  covariances  for  both  the  vertical  position  and 


velocity.  This  shows  the  type  of  tradeoff  one  expects  as  Q  varies. 


RWV(r)  Filer  portion  eetimation  error  and  coy. 


Fig.  4-2:  RWV(r)  model.  Posidon  state  estimation  errors  and  covariance  matrix  elements. 
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Fig.  4-3:  RWV(r)  model.  Velocity  state  estimation  errors  and  covariance  matrix  element. 
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Fig.  4-4:  RWV(r)  model.  Measurement  residuals  and  covariance  matrix  elements. 


4.2.2  Random  Walk  Velocity  (RWV(s))  Model  -  Range  Unavailable 

When  range  information  is  unavailable  a  spherical  coordinate  velocity  model  similar  to  Eq. 
(3-18)  is  used.  The  range  states  are  deleted  so  we  are  left  with  the  4  state  filter: 


d_ 

dt 
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(4-2) 


The  measurement  equation  is  now  a  linear  function  of  the  states  so  an  extended  Kalman  filter 
formulation  is  not  required.  The  power  spectral  density  matrix  of  the  state  process  noise  required 
for  the  discretization  equations  (2-3)  is  taken  as  (again,  see  Appendix  A) 
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The  results  of  this  "angle  only"  filter  are  shown  in  the  following  figures.  First  the  angular 
position  state  errors  and  their  estimated  covariances  are  in  Fig.  4-5,  and  the  angular  velocity  state 
errors  and  associated  estimated  covariances  are  shown  in  Fig.  4-6.  For  completeness  the 
measurement  residuals  are  shown  in  Fig.  4-7. 


RWV(i)  Filter  podtlon  ntimdion  •rrcx»  and  oov. 


Fig.  4-5:  RWV(s)  model.  Angle  state  estimation  errors  and  covariance  matrix  elements. 
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Fig.  4-6:  RWV(s)  model.  Angular  rate  state  estimation  errors  and  covariances. 
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RWV(s)  Filter  Measurement  residual*  and  cov. 
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Fig.  4-7:  RWV(s)  model.  Measurement  residuals  and  covariance  matrix  elements. 


With  the  relatively  large  noise  (1.75  mrad)  that  was  superimposed  on  the  bearing  and 
elevation  measurements,  the  filters  are  performing  about  as  well  as  could  be  expected.  Again 
the  target  maneuvers  starting  at  r=8  and  r=17  produce  brief  biases  in  the  errors  of  estimates  of 
both  the  bearing  and  bearing  rate  states.  Any  number  of  maneuver  detectors  (e.g.  Woolf  son 
(1985))  are  available  for  this  problem  but  none  are  included  in  the  simulations  to  date.  A  look 
at  the  bearing  residual  Fig.  4-7  shows  very  little  effect  of  the  maneuvers  because  of  the  high 
measurement  noise  used  in  these  examples. 


4.3  DISCUSSION 


The  simulation  results  presented  here  should  be  taken  only  as  examples  of  scenarios  that 
may  or  may  not  be  particularly  realistic.  With  the  simulation  tools  as  they  exist  now,  it  is  easy 
to  rapidly  prototype  various  filter  configurations  and  evaluate  their  performance  against  real  target 
measurement  data  from  specific  imaging  sensors.  What  is  satisfactory  for  one  type  of  sensor  may 
not  be  so  for  another.  For  example,  a  maneuver  detector  based  on  filter  residuals  may  work  well 
for  a  relatively  low  noise  sensor,  but  may  not  if  the  sensor  generates  position  data  as  noisy  as 
was  used  in  the  previous  section.  Also,  the  next  chapter  outlines  a  multi-target  algorithm  that 
uses  filter  residuals  to  determine  which  measurement  to  associate  with  which  track.  In  that  case, 
most  maneuver  detectors  based  on  residuals  will  not  work  particularly  well. 
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5.0  TRACKING  FILTER  ENHANCEMENTS 


This  chapter  outlines  a  few  enhancements  that  were  m arte  to  the  tracking  algorithms  of  the 
previous  chapter  that  may  make  them  more  applicable  to  direct-fire  support  applications.  These 
include  a  swapping  algorithm  that  allows  the  filter  to  swap  from  an  angle-only  filter  to  a  full  3 
dimensional  cartesian  coordinate  filter  when  an  active  range-finding  sensor  is  activated;  and  a 
simple  but  practical  multiple  target  tracking  algorithm  that  allows  the  fire  control  system  to 
maintain  lock  on  several  targets  in  the  field  of  view  while  another  target  is  being  engaged. 


5.1  MULTIPLE  TARGET  TRACKING 


When  a  target  is  to  be  tracked  in  a  cluttered  environment  (i.e.  one  with  potential  false 
targets)  or  when  more  than  one  target  must  be  tracked  simultaneously,  there  must  be  a  method 
in  the  algorithm  to  determine  which  track  a  particular  measurement  belongs  to,  or  if  it  is  a  false 
target  There  are  many  such  algorithms  available  in  the  literature  and  include  the  nearest 
neighbour  method,  the  probabilistic  data  association  method,  the  multiple  hypothesis  tracking 
method,  among  others,  all  of  which  are  summarized  in  Roy  (1990)  or  McMillan  and  Lim  (1990) 
or  described  in  more  detail  in  texts  such  as  Blackman  (1986)  or  Farina  and  Studer  (1985).  The 
simplest  of  these,  but  arguably  the  least  satisfactory  for  very  closely  spaced  targets,  is  the  nearest 
neighbour  method.  Nonetheless,  this  method  was  implemented  and  found  to  perform  quite  well 
given  its  simplicity.  The  method  and  some  simulation  results  are  described  in  this  chapter. 


5.1.1  Nearest  Neighbour  Method 


The  nearest  neighbour  algorithm  first  finds  all  current  tracks  that  are  reasonably  close  to 
the  measurement  by  a  process  known  as  "gating",  described  below.  If  there  are  none,  then  the 
data  is  either  a  new  target  or  a  piece  of  clutter.  In  either  case,  a  tentative  track  filter  is  initiated 
for  the  measurement  which  becomes  valid  only  if  t  lere  are  subsequent  measurements  associated 
with  it.  If,  on  the  other  hand,  there  are  one  or  more  tracks  that  are  close,  then  simply  the  track 
that  is  the  closest  is  updated  with  that  measurement. 

To  make  the  notion  of  gating  more  precise,  define  a  normalized  "distance",  d,  from  the 
location  of  the  measurement  at  time  k,  z(k),  to  the  tracking  filter’s  prediction  of  the  measure¬ 
ment,  /t(i(fclfc-l)),  for  each  track: 

d2(k\k- 1)  :=  [z(fc)  -  h(Wk- l))f  [H{Wk - 1  ))P(k\k - 1  )H  T{£(k\k - 1 ))  +  R]~l  (5.1) 

[Z(k)  -  h(£m-m 


(The  notation  here  is  the  same  as  that  of  Chapter  2.)  Note  that  this  "distance"  is  essentially  the 
length  of  the  residual  vector,  normalized  by  its  estimated  covariance.  A  gate  threshold  must  be 
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specified,  to  which  d  is  compared  to  see  if  the  measurement  is  close  to  the  filter.  This  threshold 
is  typically  around  5.  The  nearest  neighbour  algorithm  can  be  summarized  as: 

1.  Obtain  the  new  measurement,  z(k). 

2.  Propagate  all  current  tracks  to  the  new  time  (i.e.,  compute  £(k\k-\),  P(k|jk-1)). 

3.  Compute  d  as  in  (5-1)  for  each  track. 

4.  If  d  of  all  tracks  >  gate  threshold  then  initiate  a  new  tentative  track  filter. 

5.  Otherwise  update  the  track  having  the  minimum  d  with  the  measurement 

Tracks  can  be  deleted  if  they  are  not  updated  with  a  new  sensor  measurement  after  a  certain 
period  of  time  or  a  certain  number  of  scans. 


5.2  FILTER  SWAPPING  ALGORITHMS 


Typically  the  fire  control  system  searches  for  targets  in  a  passive  mode  with  an  optical  or 
infrared  imaging  sensor  and  tracks  these  targets  with  angle-only  filters.  When  a  target  is  to  be 
engaged,  an  active  ranging  sensor  is  turned  on  and  information  from  the  third  spatial  dimension 
becomes  available  to  the  tracking  filters.  It  would  be  preferred,  then,  to  use  one  of  the  cartesian 
coordinate  target  models  since  target  motions  are  more  naturally  described  in  cartesian 
coordinates.  If  for  some  reason  the  target  is  disengaged  but  still  is  to  be  tracked,  the  filters  will 
have  to  revert  to  a  passive,  angle-only  tracking  mode.  This  section  describes  a  method  that  can 
be  used  to  swap  a  2  (spatial)  dimensional  passive,  angle-only  tracking  filter  to  a  3  dimensional 
active,  cartesian  one  when  range  information  becomes  available  and  vice  versa  when  it  ceases. 


5.2.1  Passive  to  Active  Filter  Swap 


In  essence,  the  filter  swapping  technique  uses  the  current  state  estimate  and  covariance 
matrix  to  initialize  another  Kalman  filter  that  uses  the  new  set  of  dynamical  equations.  To  swap 

from  an  angular  filter  (with  current  state  vector  [bib  t  b  e  ]r)  to  a  rectangular  one  (with 

a  new  state  vector  [  pXo  Pyg  Pzg  ^xg  ^yg  ^zg  ^xg  &yg  &zc  lr)  when  a  valid  piece  of  range, 

range  rate  or  range  acceleration  (zr ,  zt,  or  zr)  data  becomes  available,  the  following  steps  can 
be  initiated: 
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1.  Compute  the  Line-of-sight  to  Geographic  DCM,  CG  ,  using  the  current  best  estimate 
of  bearing  and  elevation  available  from  the  filter: 


cos  2  cos S 
cosi  sin  b 


-sin£ 


-sin£  sin£  cos B 
cos  6  sin  #  sin  S 
0  cos  S 


2.  Estimate  the  geographic  relative  position  vector  from  the  measured  range  and  filtered 
bearing  and  elevation: 
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3.  Estimate  the  line-of-sight  relative  velocities  from  the  estimated  angular  rates  and  the 
measured  range  and  range  rates  (if  range  rate  not  available,  use  zt  =  0): 


^  := 


1L) 


XL) 


Zr  b  cos  2 

A 

~zr  k 


38 


4.  Convert  the  line-of-sight  velocities  to  geographic  velocities: 


5.  Estimate  the  line-of-sight  accelerations  from  the  estimated  angular  accelerations  and 
the  measured  range  and  range  rates  (if  range  acceleration  not  available,  use  zf  =  0): 
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6.  Convert  the  line-of-sight  accelerations  to  geographic  accelerations: 


The  corresponding  swapping  of  exact  covariance  matrix  information  would  be  very 
complex  (due  to  the  nonlinearity  of  the  coordinate  transformations)  and  of  minimal  value  since 
the  covariance  matrix  quickly  converges  to  its  steady  state  values.  The  following  is  an 
approximation  that  will  initialize  a  new  diagonal  covariance  matrix  from  the  diagonals  of  the  old 
and  which  should  result  in  relatively  modest  covariance  element  transients  following  a  filter 
swap: 
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1.  Extract  the  square  roots  of  the  diagonals  of  the  existing  covariance  matrix: 
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It  should  be  noted  that  a  passive  filter  may  not  have  range  states  and  consequently  no 
corresponding  covariance  elements  (Prr  ,  Ptf.,  Pff).  The  filter  obtains  estimates  for  these 
parameters  from  a  configuration  file. 

2.  Approximate  the  corresponding  L  frame  standard  deviations  via: 


°Pxl 

[L) 

<*r 

For  position  states: 

°Pyl 

**Pzl 

SO 

zr  ab  |cos2| 

.  Zr°‘  . 

ov 

VXL 

[L) 

For  velocity  states: 

Ov 

VYL 

Ov 

VZ L 

SO 

zr  ob  |cos2| 

zr*t 

For  acceleration  states: 
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(Compare  these  transformations  with  Eqn.  (3-8)  which  related  the  bearing,  elevation  and 
range  rates  to  the  velocites  in  L  coordinates.  Such  a  transformation  takes  the  components 
of  the  covariance  matrix  representing  the  uncertainty  in,  say,  the  angular  position  estimate 
(or,  Ojj,  ce)  and  uses  the  measured  range  and  elevation  to  convert  these  to  rectangular 
coordinates.  The  uncertainty  in  range,  <5^  transforms  directly  to  GpxL*  since  the  axis 
points  along  the  range.  The  uncertainty  in  elevation,  cg,  multiplied  by  the  range  gives  the 
uncertainty  in  the  vertical  cartesian  direction  ZL.  Similarly  for  the  bearing  and  the  YL 
direction.  The  same  transformation  is  applied  to  the  angular  velocity  and  acceleration 
uncertainties.) 

3.  Convert  these  to  G  frame  coordinates  with  the  current  DCM: 
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For  velocity  states: 
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For  acceleration  states: 

^aYG 
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rv 

■"a 

Gfli Z 

aTJG 
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4.  Initiate  the  covariance  matrix  as: 

/>(G)  :=  diag  [QPx(}'  apYG  ♦  °Pzc  ’  °vxc  ’  °vrc  ’  °vzg  ’  °aXG  ’  °aYG  ’  °azc  ^ 


5.2.2  Active  to  Passive  Filter  Swap 

When  the  system  loses  ranging  information  to  the  target,  the  filter  will  coast  through  a  few 
updates  using  the  predicted  range  measurement  as  new  data.  After  a  certain  number  of 
successive  updates  with  no  new  range  data,  the  filter  will  revert  to  a  bearing  and  elevation  angle 
tracker.  The  initial  conditions  for  the  new  filter  are  derived  from  the  last  state  and  covariance 
matrix  estimates  tom  the  cartesian  filter  in  a  manner  similar  to  that  above.  To  obtain  the  new 

state  vector  [  b  e  b  e  b  e  ]T  from  the  old  [  pXc  PyG  PzG  ^XG  ^Yc  ^ZG  ^XG  ^Yc  ^ZG 
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1.  Compute  the  Geographic  to  Line-of-sight  DCM,  Cq  ,  using  the  current  best  estimate 
of  bearing  and  elevation  available  from  the  filter: 

f  :=  fac  +  Pyg  + 

b  :=  tan'1  WyJHJ 
i  :=  sin'1  {-PzJt) 

cos  &  cos  b  cosS  sin 6  -sinS 
-sinh  cos£  0 
sin£  cosh  sin£  sine  cost?  _ 

2.  Convert  the  geographic  velocities  to  line-of-sight  velocities: 
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3.  Estimate  the  angular  rates  from  the  line-of-sight  velocities: 
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4.  Convert  the  geographic  accelerations  to  line-of-sight  accelerations: 
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5.  Estimate  the  angular  accelerations  from  the  line-of-sight  accelerations: 


A 

? 

&xL  +  +  n2cos2s 

5 

ss 

dY  /(f cos S)  -  2(f/P)b  +  2  b  £ tan£ 

L$ 

A 

i 

-&z,W  ~  2(^)  e  ~  b2cosismi 

Id 

The  covariance  matrix  for  the  new  filter  can  be  initiated  from: 

1.  Extract  the  square  roots  of  the  diagonals  of  the  existing  covariance  matrix: 


p 

GPXC 

^PXG  PXG 

For  position  states: 

gPyg 

:= 

Pyg  Pyg 
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For  velocity  states: 
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:= 
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For  acceleration  states: 

®aYG 

:= 

\j^aYG  aYG 

8 
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2.  Convert  these  to  L  frame  coordinates  with  the  current  DCM: 


For  position  states: 


For  velocity  states: 
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For  acceleration  states: 
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3.  Approximate  the  corresponding  spherical  frame  standard  deviations  via: 
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apYL/(fcosi) 

For  position  states: 
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4.  Initiate  the  covariance  matrix  as: 


P  :=  diag  [ab2,  a2,  ab2,  a*2,  ...] 


5.3  SIMULATIONS  OF  ENHANCED  TRACKER  ALGORITHMS 


This  section  shows  some  numerical  results  of  the  filter  swapping  algorithms  described  in 
the  previous  sections.  Consider  the  trajectories  of  the  two  targets  shown  in  Fig.  5-1.  Thi< 
represents  two  targets  about  20  m  apart  and  about  1600  m  from  the  sensor  moving  at  54.0  km/hr. 
Initially  they  are  travelling  away  from  the  sensor  but  then  turn  and  head  towards  the  sensor. 
Assume  passive  sensors  are  used  except  for  the  portions  of  the  trajectories  shown  in  the 
rectangles  (from  r=15  to  f=20  for  Target  1  and  from  r=17  to  r=19  for  Target  2)  when  a  ranging 
sensor  is  turned  on. 


5.3.1  Simulations  of  Filter  Swapping 


To  demonstrate  the  performance  of  the  filter  swapping  algorithms,  consider  only  a  single 
target  (Target  1  of  Fig.  5-1,  the  same  target  in  example  of  the  previous  chapter)  which  is  tracked 
with  a  4  state  bearing-elevation  angle  filter  for  15  seconds.  At  t=15  seconds,  valid  range 
measurements  are  available  for  the  next  5  seconds  and  then  cease  for  the  final  5  seconds  of  the 
engagement.  The  Kalman  filter  matrices  and  parameters  were  the  same  as  those  used  in  the 
previous  chapter.  Shown  in  Fig.  5-2  below  is  the  first  element  of  the  state  vector.  For  the  first 
15  seconds  and  the  final  5,  it  represents  the  filter  estimate  of  the  target  bearing  angle  in  radians. 
Between  r=15  and  f=20,  it  is  the  estimate  of  pXQ,  the  northerly  distance  to  the  target  in  meters. 
The  truth  data  for  these  are  shown  in  Fig.  5-3.  Since  the  scales  of  the  different  portions  of  the 
graph  are  not  compatible,  the  next  two  figures,  Fig.  5-4  and  Fig.  5-5,  show  the  distinct  regions 
of  interest.  The  corresponding  covariance  estimates  are  also  shown  in  Fig.  5-6  and  Fig.  5-7. 
Again  since  the  scales  of  the  different  segments  of  the  plots  are  not  compatible,  two  figures  are 
used. 
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Fig.  5-1:  True  x-y  trajectories  of  two  targets.  Ranging  information  is  assumed  to  be  available 
in  the  portions  of  the  trajectories  enclosed  by  rectangles. 


Fig.  5-6:  Covariance  of  first  state  Fig.  5-7:  Covariance  of  first  state  (XG 

(bearing)  when  range  is  not  position)  when  range  is 

available.  Notice  small  available  from  f=15  to  f=20. 

transient  after  swap  at  r=20. 


The  covariance  of  the  pXG  state  in  Fig.  5-7  shows  some  initial  transients  after  it  was  formed 
when  range  became  available  at  r=15.  When  range  information  ceased  at  f=20,  the  cartesian  filter 
tried  to  coast  along  for  the  next  three  updates  (0.3s)  with  no  range  measurements  and 
consequently  the  covariance  estimates  started  to  diverge  rapidly  at  r=20.  The  filter  was  then 
swapped  back  to  the  bearing-elevation  filter  and  the  corresponding  covariance  element  became 
that  of  the  bearing  state  again,  Fig.  5-6,  after  some  brief  transients  after  t=20. 

5.3.2  Simulations  of  Multi-Target  Nearest  Neighbour  Tracking  With  Filter  Swapping 

The  simulation  system  previously  described  was  modified  to  include  the  nearest  neighbour 
multitarget  algorithm.  Sensor  data  files  of  closely  spaced  targets  were  generated  and  the  files 
were  merged  and  sorted  by  time  stamp  with  no  indication  as  to  which  target  generated  which 
return.  The  algorithm  associates  a  new  measurement  to  the  closest  existing  track  (if  there  is  one 
within  the  gate  threshold)  or  it  spawns  a  new  filter  for  a  measurement  lying  outside  all 
thresholds.  When  the  filters  were  properly  tuned,  the  algorithm  would  eventually  settle  on  one 
filter  for  each  target  with  only  the  occasional  misassociation.  Clutter  points  generated  at  random 
would  also  spawn  new  filters,  but  without  subsequent  measurements  being  associated  with  them, 
these  filters  would  quickly  die. 
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Fig.  5-8  below  shows  the  true  and  measured  bearings  to  the  two  targets  of  Fig.  5-1. 
Likewise  Fig.  5-9  shows  the  true  and  measured  elevations.  The  true  and  measured  ranges  shown 
in  Fig.  5-10  require  a  bit  of  explanation:  the  measured  range  is  negative  in  the  portions  of  the 
trajectory  when  the  ranging  sensor  is  turned  off.  That  is  how  the  simulator  distinguishes  between 
a  passive  and  an  active  sensor.  The  performance  of  the  nearest  neighbour  algorithm  is 
represented  in  the  top  chart  of  Fig.  5-11.  A  spike  in  the  plot  for  Filter  1,  for  example,  means  that 
it  was  incorrectly  updated  with  a  measurement  from  target  2.  Fig.  5-12  shows  the  second  state 
variable  from  Filter  1  (which  is  the  elevation  state  when  it  was  a  bearing-elevation  angle  filter, 
and  is  the  east  position  state  when  it  is  swapped  to  a  cartesian  filter  at  f=15).  The  bottom  portion 
of  the  same  figure  shows  the  corresponding  state  error  and  its  covariance.  The  Kalman  filter 
matrices  and  parameters  remained  same  as  in  the  previous  example,  and  the  "gate  thresholds" 
were  chosen  as  6.0  for  a  cartesian  filter  and  4.5  for  a  bearing-elevation  filter.  Filter  1  is  spawned 
immediately  and  takes  all  the  measurements  from  both  targets  and  associates  them  with  the  same 
filter  for  the  first  2  seconds  or  so.  Eventually  the  filter  converges  sufficiently  to  recognize  there 
are  two  distinct  targets  and  spawns  a  new  filter  for  Target  2  at  about  t=l.  From  then  on,  the  two 
targets  are  distinguished  most  of  the  time,  with  occasional  misassociations  occurring  at  the  times 
of  the  spikes.  The  lower  graph  in  Fig.  5-11  represents  the  times  that  the  filters  were  swapped 
as  the  range  information  became  available  or  ceased:  Filter  1  was  swapped  from  a  BE  (bearing- 
elevation)  filter  to  an  XYZ  (cartesian)  filter  when  the  range  information  of  Target  1  was  valid 
from  r=15  and  reverted  to  the  BE  filter  shortly  after  r=20.  Likewise  for  Filter  2  when  the  range 
of  Target  2  was  available  from  r=17  to  r=19. 


Measured  Bearing  Angles,  Targets  1  &  2 
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Fig.  5-8  Measured  and  true  bearing 
angles,  two  closely  spaced 
targets 


Measured  Elevation  Angles,  Targets  1  &  2 


Fig.  5-9  Measured  and  true  elevation 
angles,  two  targets 
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Measured  Range,  Targets  1  &2 
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Fig.  5-10  Measured  and  true  range,  two 
targets.  (Negative  range 
indicates  it  is  unavailable  to 
the  filter.) 


Fig.  5-11  Times  of  misassociations  (top 
figure)  and  times  of  filter 
swaps  (below). 


Filter  1 ,  state  2  estimate,  error  and  covariance 


Fig.  5-12:  Elevation  state  estimate  from  Filter  1,  showing  misassociations  &  swap 


6.0  KALMAN  FILTERS  FOR  MUZZLE  REFERENCE  SYSTEMS 


6.1  INTRODUCTION 


This  chapter  describes  the  modelling  and  simulation  work  that  was  done  with  the  goal  of 
applying  Kalman  filtering  to  a  dynamic  muzzle  referencing  system.  The  DMRS  hardware 
consists  of  a  small  laser  transceiver  mounted  on  or  near  the  mantlet  of  the  gun.  A  small  mirror 
is  mounted  on  the  muzzle  to  reflect  the  laser  beam  back  to  the  receiver.  As  the  barrel  flexes  due 
the  motion  of  the  vehicle  over  rough  ground,  the  angle  of  the  laser  return  will  deviate 
proportionally.  This  angle  can  be  detected,  sampled  and  made  available  to  the  fire  control 
system.  Signal  processing  is  then  required  to  process  these  angles  and  predict  when  the  barrel 
flex  will  be  close  enough  to  zero  so  that  the  shell  can  be  launched  through  a  barrel  that  is 
relatively  straight.  A  preliminary  Kalman  filter  design  is  presented  in  this  chapter  and  is  shown 
to  perform  this  prediction  quite  well. 

First,  the  gun  barrel  dynamics  will  be  modelled  by  considering  the  first  four  modes  of 
vibration  of  a  long,  flexible,  hollow  tube  with  the  physical  dimensions  of  a  tank  gun.  Then,  these 
equations  are  transformed  to  a  series  of  state  equations  that  can  be  used  for  simulations.  The 
power  spectrum  of  the  simulated  gun  is  compared  with  that  obtained  from  actual  field  data. 
Finally,  a  Kalman  filter  is  designed  with  a  subset  of  these  state  equations  and  its  ability  to  predict 
the  muzzle  pointing  angle  is  tested  on  actual  data. 

There  has  been  some  ground  work  that  used  Kalman  filters  in  dynamic  muzzle  referencing 
systems.  In  particular.  Levin  (1978)  used  a  simple  second  order  model  for  an  analog  Kalman 
filter  and  Baran  et  al.  (1987)  described  an  algorithm  that  first  estimated  the  gun  barrel  model  (by 
finding  the  coefficients  of  an  auto-regressive,  moving  average  process)  with  the  first  few  seconds 
of  barrel  dynamics  measurements  and  then  used  these  coefficients  in  a  Kalman  filter  to  predict 
the  muzzle  angle  20  ms  later.  The  encouraging  results  led  them  to  demonstrate  the  system  in 
real  time  on  an  M1A1  tank  ,  as  described  in  Brosseau  et  al.  (1990).  The  methods  described  in 
this  chapter  offer  another  alternative  to  developing  the  gun  barrel  model  for  the  predictive 
Kalman  filter;  that  is,  starting  from  the  modal  equations  of  a  long,  uniform,  hollow  tube  that 
approximates  the  gun  barrel. 


6.2  TRANSVERSE  VIBRATIONS  OF  LONG  HOLLOW  BEAMS 


The  theory  of  vibration  of  uniform  beams  is  a  well-studied  area.  In  general,  however  a  gun 
barrel  is  not  uniform.  It  has  several  different  cross  sections  along  its  length  and  is  very  often 
tapered  at  some  points.  Closed  form  solutions  for  the  modal  shapes  and  frequencies  are  not 
known  in  general  and  require  extensive  modelling  and  numerical  simulation  efforts  to 
approximate.  Such  modelling  work  is  done  by  gun  designers  and  is  presented  in  detail  by  Gast 
(1987),  for  example. 

For  the  purposes  of  this  report,  however,  we  will  assume  the  gun  barrel  can  be  modelled 
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as  a  uniform  hollow  tube  that  is  free  at  the  muzzle  end  and  is  either  fixed  or  hinged  (depending 
on  the  operation  of  the  gun  control  system)  at  the  mantlet  end.  The  development  in  this  chapter 
uses  standard  theory  of  free,  undamped,  transverse  vibrations  of  linearly  elastic  beams,  also 
known  as  Bemoulli-Euler  beam  theory  (e.g.  Craig  (1991)).  There  are  a  number  of  assumptions 
in  Bemoulli-Euler  theory: 

a)  the  longitudinal  axis  of  the  beam  undergoes  no  extension  or  contraction; 

b)  transverse  shear  deformation  is  neglected; 

c)  the  material  is  linearly  elastic; 

d)  the  beam  has  uniform  cross-section;  and 

e)  the  slope  of  the  beam  remains  small. 

These  assumptions  are  not  unreasonable  for  an  approximation  of  the  gun  barrel  that  is  suitable 
for  a  short  term  predictive  model  for  the  Kalman  filter. 

Denote  the  displacement  of  the  centerline  of  the  beam  from  its  at-rest  centerline  as  l }(l,t). 
Note  that  T|  is  a  function  of  both  the  distance  along  the  beam,  /,  and  time,  t.  See  Fig.  6-1  for  a 
pictorial  representation  of  the  beam.  Using  standard  Newtonian  mechanics  and  the  above 
assumptions,  it  can  be  shown  that  t]  satisfies 


/£/32t1(/,/)A 


V 


a/2 


+  pA 


&r\(M 

dt2 


(6-1) 


where  f(l,t)  is  an  external  applied  force,  E  is  the  elastic  modulus  of  the  material,  p  is  its  density, 
A  is  the  cross-sectional  area  and  7  is  the  transverse  moment  of  inertia  of  the  beam. 


Fig.  6-1:  Schematic  of  generic  flexible  beam 
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There  are  different  boundary  conditions  for  this  differential  equation  depending  on  the 
physical  constraints  placed  on  the  beam  for  fixed,  free  or  hinged  ends: 


f  nd'J)  =■  o 


Fixed  end: 


(6-2) 


Free  end: 


a/2 

Hi* 


=  o 


J  Ul. 


=  0 


Jw, 


(6-3) 


r\de.t)  =  0 


Hinged  end: 


=  0 


(6-4) 


Specifying  the  type  of  support  at  each  end  of  the  beam  will  provide  4  boundary  conditions.  To 
obtain  a  solution  for  the  unforced  (f[l,t)=0 )  part  of  (6-1),  we  must  make  an  assumption  on  the 
form  of  the  solution  T It  is  assumed  that  the  temporal  and  spatial  portions  of  rj  can  be 
separated  as  the  product  of  a  modal  shape  function,  0(7),  (that  depends  only  on  the  distance  0 
and  a  simple  harmonic  function,  h(t)  =  c0  cos(cor  -  a),  that  depends  only  on  time: 

rj  (/,r)  =0(0  h(t)  =  c0  0(0  cos(co  t  -  a)  (6-5) 


where  c0  is  a  constant  From  this,  we  see  that 


=  -c0  co2  0(0  cos(o)  t  -  a) 
dt2 


and 
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cQ  cos(co/  -  a) 


a4r|(/,r)  =  Affl 

a/4  diA 

Substituting  these  in  (6-1)  and  setting  f(l,t)=Q  yields 

os(o)f-a)  -  p  A  c0g)2$ (/)cos(o)f-a)  =  0 
dl 4 


or  by  defining 


X4  :=  PA(°2  (6-6) 

El 

then 

-  X4<j)(/)  =  0  (6-7) 

dl4 

A  general  solution  to  (6-7)  is  well  known  to  be 

<)>(/)=  c1sinh(X/)  +  c2  cosh(X/)  +  c3sin(X/)  +  c4  cosQJ)  (6-8) 

where  we  have  5  unknowns  (the  4  amplitude  coefficients  and  the  eigenvalue  X).  The  next  two 
subsections  use  different  sets  of  boundary  conditions  for  two  configurations  of  the  beam:  the 
hinged-free  beam  and  the  fixed-free  beam. 


6.2.1  Hinged-Free  Beam:  Modal  Shapes  and  Frequencies 


Assume  the  beam  is  hinged  at  /  =  0  and  free  at  l  =  L.  Now  apply  the  4  known  boundary 
conditions.  For  a  beam  hinged  at  /  =  0  boundary  conditions  (6-4)  imply 


d26 
dl 2 


6(0)  =0  =>  c2  +  c4  =  0 
=  0  =>  X2(c2  -  =  0 

J/«o 


(6-9) 
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and  the  free  boundary  conditions  at  /  =  L  imply,  from  (6-3), 


=  0  =>  X2  (Cj  sinhXL  +  c2  coshXL  -  c3  sinXL  -  c4  cosXL)  =  0 
/*L 

(6-10) 

=  0  =>  X3  (c^ coshXL  +  c2sinhXL  -  c3 cosXL  +  c4  sinXL)  =  0 
/«L 

To  get  a  nontrivial  (X*0)  solution  to  (6-9)  we  see  that 

c2  =  c4  ~  0 

Substituting  this  result  in  (6-10)  will  result  in 

X2(c,  sinhXL  -  c3  sinXL)  *  0 

(6-11) 

X3  (Cj  coshXL  -  c3  cosXL)  =  0 
which  can  be  written  in  matrix  form  as: 


X2  sinhXL  -X2  sinXL 

'cl 

0 

X3  coshXL  -X3  cosXL 

53. 

0 

Again  for  this  to  have  a  non-trivial  solution,  we  can  set  the  determinant  of  this  matrix  to  0  and 
solve  for  XL: 

X2  sinhXL  -X2  sinXL 

=  0 

X3  coshXL  -X3  cosXL 
=>  sinhXL  cosXL  -  coshXL  sinXL  =  0 


d2( fr 
dl2 

dl3 
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This  equation  has  many  solutions,  a  few  of  which  were  found  numerically  and  are  listed  here: 


XjL  =  3.9266 
k^L  =  7.0686 
X3 L  =  10.2102 
XaL  =  13.3518 


Now  from  the  definition  of  X  in  (6-6),  we  can  get  the  modal  frequencies 


(0 


2 


x4l4ei 

pAL4 


=> 


0 )t-  =  (kpy 


El 


p  AL4 


(6-12) 


(6-13) 


Table  6-1:  Approximate  physical  parameters  of  a  105mm  gun  barrel 

L 

=  Length  of  barrel  =  5  m 

dQ 

=  Outer  diameter  =  0.125  m 

di 

=  Inner  diameter  =  0.105  m 

A 

=  Cross-sectional  area  =  (n/4 )(d02  -  d2)  *  0.00361  m2 

I 

=  Cross-sectional  moment  of  inertia  =  (n/64)(d4  -  d4)  =  6.02x1  O'6  m4 

E 

=  Elastic  modulus  of  steel  =  2.068X1011  N/m2 

P 

=  Density  of  steel  =  8000  kg/m3 

For  a  generic  105mm  tank  gun  barrel,  the  approximate  physical  parameters  in  Table  6-1 
can  be  used.  Combining  these  numerical  values  with  the  values  of  kjL  from  (6-12)  and 
substituting  in  (6-13)  yields  the  modal  frequencies  for  the  hinged- free  hollow  gun  tube: 


©!  =  (3.92662)  (8.3)  =  128  rad/s  =  20.4  Hz 
co2  =  (7.06862)  (8.3)  =  415  rad/s  =  66  Hz 
to3  =  (10.21022)  (8.3)  =  865  rad/s  =  138  Hz 
©4  =  (13.35182)  (8.3)  =  1480  rad/s  =  235  Hz 


To  obtain  the  corresponding  modal  shape  functions,  we  return  to  finding  the  amplitude 
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coefficients  cv  From  the  first  equation  of  (6-11)  we  can  express 

sinAZ, 

Cl  sinhAX.  3 


and  since  we  previously  established  that  c2  =  c4  =  0  then  we  get  from  (6-8)  the  modal  shape 
functions  for  the  hinged-free  tube: 


♦|d>  =  c3 


sinX^Ll.  + 


^  sin X^  ^ 


sinhXX 
V  *7 


sinhX^J. 


These  modal  shape  functions  (with  c3  arbitrarily  set  to  1)  are  plotted  in  Fig.  6-2. 


Modal  shapes:  Hinged-Free 


(6-15) 


Fig.  6-2:  Modal  Shape  Functions  of  Hinged-Free  Hollow  Tube 


6.2.2  Fixed-Free  Beam:  Modal  Shapes  and  Frequencies 

In  the  second  situation,  we  assume  the  beam  is  fixed  at  /  =  0  and  free  at  l  =  L.  For  this 
case  the  boundary  conditions  at  /  =  0  come  from  (6-2): 
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<J>(0)  =0  =>  c2  +  c4  -  0 

=  0  =>  X.(Cj  +  c3)  =  0 
ImO 


and,  as  before,  the  free  boundary  conditions  at  /  =  L  imply  from  (6*3) 


d2( j> 
dl 2 


■  0  =>  X2  (Cj  sinhAX  +  c2coshA£  -  c3sinXL  -  c4cosAjL)  =  0 


J  l*L 


d\ > 
dl 2 


=  0  =>  X3  (^coshXL  +  c2sinhAL  -  c3cosXL  +  c4sinXL)  =  0 


J/*L 


These  4  equations  can  be  written  in  matrix  form  as: 


0  1  0  1  1 

ci 

0 

o 

€< 

© 

<< 

'2 

0 

X2  sinhAZ,  X2  coshAZ  -X2  sinAL  -A.2  cosAL 

c3 

0 

A,3  coshAL  X3  sinhXL  -X3  cosXL  -X3  sinAL 

0 

(6-16) 


As  before,  a  nontrivial  set  of  solutions  to  this  set  of  equations  is  obtained  if  the  determinant  of 
the  matrix  is  set  to  zero: 


0  10  1 
X  0  X  0 


X2sinhAJL  X2coshAL  -X2  sinXL  -A,2  cos  XL 
A.3  coshAL  A.3  sinhAX  -A3cosAL  -A.3  sinAL 


=  0 


=>  cosAL  coshAL  +  1=0 

Numerical  solution  of  this  equation  leads  us  again  to  several  solutions,  one  for  each  mode: 

XxL  =  1.8751 
X^L  =  4.6941 
X^L  =  7.8548 
XaL  =  10.996 
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so  substituting  these  (6-13)  yields  the  modal  frequencies  for  the  fixed-free  hollow  gun  tube: 


=  (1.8752)  (8.3)  =  29  rad/s  =  4.6  Hz 
co2  =  (4.6942)  (8.3)  =  183  rad/s  =  29  Hz 
co3  =  (7.8552)  (8.3)  =  512  rad/s  =  82  Hz 
co 4  =  (10.9962)  (8.3)  =  1004  rad/s  =  160  Hz 


(6-17) 


To  obtain  the  corresponding  modal  shape  functions  of  the  fixed-free  beam,  we  need  the 
amplitude  coefficients  cf.  It  is  straightforward  to  show  from  (6-16)  that 

c4  =  -c2 

c3  =  ~c\ 


c\  ~  ~c2 


coshX^L  +  cos  XjL 
sinhA/.  +  sinXjL 


so  that  by  leaving  c2  arbitrary  and  substituting  these  coefficients  in  the  general  modal  shape 
equation  (6-8),  then  the  mode  shapes  for  the  fixed-free  beam  are  (and  as  plotted  in  Fig.  6-3): 


$j(  0  =  c2 


coshA,jL.i-  -  cosXjL—  -  (sinhAjL-L  -  sinX^L-^l 


where  ri 


coshX^  +  coshXjL 
sinhXjL  +  sinXjL 


(6-18) 


6.3  FORMULATION  OF  STATE  EQUATIONS 


With  the  assumed  form  of  the  solution  of  the  main  differential  equation  (6-1)  describing 
the  dynamics  of  a  transversely  vibrating  beam  having  been  established  in  the  previous  section 
as  a  product  of  the  modal  shape  functions,  (6-15)  or  (6-18),  with  time-varying  simple  harmonic 
functions,  it  is  possible  to  derive  a  set  of  state  equations  that  describe  the  relative  contributions 
of  these  mode  shapes  at  any  point  along  the  beam,  at  any  time  and  given  any  forcing  function 
and/or  initial  conditions.  These  state  equations  then  form  the  basis  of  the  prediction  filter. 


6.3.1  Development  of  the  State  Equations 


The  method  described  in  Sasiadek  and  Srinivasan  (1989)  was  used  to  derive  the  state 
equations  and  is  summarized  here.  The  beam  is  considered  to  be  of  the  fixed-free  configuration 
of  Section  6.2.2  so  the  modal  shape  functions  of  (6-18)  are  used,  but  an  additional  degree  of 
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2.0 


Modal  shapes:  Fixed-Free 


1.0 


£  0.0 

I 

§ 

£ 

•o 

-1.0 


-2.0 

0.0  1.0  2.0  3.0  4.0  5.0 

Length  (m) 

Fig.  6-3:  Modal  Shape  Functions  of  Fixed-Free  Hollow  Tube 

freedom  is  introduced  at  the  "fixed"  end  by  allowing  that  end  to  rotate  by  an  angle  0  (see 
Fig.  6-1).  In  this  way  we  could  introduce  a  control  action  at  the  fixed  end  to  drive  the  beam  to 
any  desired  pointing  angle.  In  essence  this  is  a  reasonable  representation  of  the  gun  barrel 
attached  to  the  gun  control  system. 

The  assumed  mode  solution  describes  the  deflection  of  any  point  of  the  beam  as  the 
product  of  the  modal  shape  functions,  §t(l)  from  (6-18),  with  arbitrary  temporal  functions, 
which  in  vector  form  is  denoted  by 

ri(/,r)  =  e(tf  m  =  <t>,(0)  (&19) 

i-1 

If  we  choose  our  state  vector  xc  as 

xc(t)  :=  j  6(t)  ex(t)  e2(t)  -  en(t)  ]7 

then  Sasiadek  and  Srimvasan  (1989)  show  that  the  dynamic  equation  representing  the  undamped 
transverse  vibrations  of  the  beam  pictured  in  Fig.  6-1  ’s 
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MXc+Kxc=f 


(6-20) 


where 


M  := 


mQ  mt 


w 


o 


? 


and  where 


mQ  :=  pALrl'i 


Li 

we<j>  :=  P AfoiDdl 
o 

L 

:=  pAfe(I)$T([)dl 


K,  :=Ei\im  *mTdi 

*  i  a/2  a/2 

x  :=  applied  torque  at  the  hinged  end 

Note  that  the  second  derivatives  of  the  modal  shape  functions  are  obtained  from  (6-18)  as 


a26;(/) 


coshXjL-i  +  cosXjL-i  -  Tj  sinhX^L-i  +  smkJLJ-  j 
L  L  ^  L  L  J_ 


Since,  in  reality,  the  oscillations  are  not  undamped,  we  introduce  a  damping  matrix  into  (6-20) 
denoted  by  D  :=  diag  [  dQ  dl  d2  ...  dn  ]  so  that  the  damped  dynamic  equation  becomes 

M  X  +  D  x  +  K  xc  -  f 


or  in  more  convenient  state  space  form: 
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>c 

-M~lD 

-rrlK 

*c 

A/"1 

+ 

*c_ 

m 

[°] 

xc_ 

J01  _ 

(6-21) 


where  [/]  represents  the  (n+1)  by  (n+1)  identity  matrix  and  [0]  is  the  similarly  dimensioned  zero 
matrix.  The  gun  barrel  dynamics  as  specified  in  this  last  equation  are  now  of  the  form 
x  =  Fx  +  Gu  that  can  readily  be  modelled  in  a  Kalman  filter. 


6.3 2  Numerical  Calculations  and  Simulations 


To  make  the  state  equation,  (6-21),  more  concrete,  we  can  compute  some  numerical  values. 
Commercial  numerical  software  was  used  to  compute  the  various  integrals  for  terms  of  the  state 
dynamics  matrices  and  the  inverse  and  products  of  the  various  matrices.  The  physical  gun 
parameters  of  Table  6-1  were  used  to  calculate  the  first  4  modes.  The  integrals  evaluated  to 

mQ  =  1204.28 


ntc 


410.98 

144.5 

0 

0 

0 

65.62 

it 

af 

0 

144.5 

0 

0 

23.47 

0 

0 

144.5 

0 

11.89 

.  0 

0 

0 

144.5_ 

K*  = 


123054 

-18.22 

42.66 

-47.93 

-18.22 

4.833xl06 

-1114.5 

1223.54 

42.66 

-1114.5 

3.79xl07 

-2230.7 

-47.93 

1223.54 

-2230.7 

1.455xl0f 

and  when  the  damping  matrix  was  chosen  as 

D  =  diag  [2000  0.1  0.1  0.1  O.l] 

then  the  continuous  time  state  equations,  (6-21),  with  an  abbreviation  of  notation,  are 

x  =  Fx  +  Gf 


x  := 


-I* 


e1  &2  0  £-7  e. 


1  e2  e3  e4 
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where  the  coefficient  matrices  are 


F  matrix: 


1 . 665e+07 
4 . ?3?e+07 
7 . 563e+06 
2. 7Q4e+06 
2.377e+06 
0 
0 
0 
0 
0 


C  matrix: 


1.391 

-3.95651 

-0.631614 

-0.225851 

-0.114447 

-3.95651 

11.2607 

1.79654 

0.642401 

0.325529 

-0.631614 

1.79654 

0.293719 

0.102552 

5.196719e-02 

-0.225851 

0.642401 

0.102552 

4 . 358950e-02 

1 . 858225e-02 

-0.114447 

0.325529 

5.196719e-02 

1 . 858225e-02 

1 . 633524e-02 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-2782  0.39  5  6.316e-02  2.258e-02  1.144e-02 

7913  -1.12  -0.179  -6. 424e-02  -3.255e-02 

1263.2  -0.179  -2 . 93e-02  -1.02Se-02  -5.196e-03 

451.7  -6. 424e-02  -1.025e-02  -4.358e-03  -1.858e-03 
228.9  -3.255e-02  -5.196e-03  -1.858e-03  -1.633e-03 


0  4 . 87e+05  3.052e+06  8.558e+06 
0  -1 . 38e+06  -8.682e+06  -2.434e+07 
0  -2.21e+05  -1 . 41 9e+06  -3.886e+06 


0 

0 

0 

0 

0 

0 

0 


-79049  -4 . 956e+05  -1.651e+06 
-40057  -2 . 511e+05  -7.041e+05 


The  time-varying  magnitudes  of  the  various  states  are  dependent  on  the  initial  conditions 
and  input  forcing  function  /.  Computer  simulations  were  conducted  on  the  above  barrel  model 
with  a  relatively  large  (2  cm)  initial  amplitude  on  mode  1  (state  e{)  and  no  input  forces,  /= 0. 
The  output,  y(t),  was  chosen  as  the  displacement  at  the  free  end  of  the  beam: 

4 

yd)  :=  n<M  «  £(«,</)*,<«) 

i=l 

=  2 ex{f)  -  le2{t)  +  2 e3(t)  -  2e4(r) 

=  [  0  0  0  0  0  0  2  -2  2  -2  ]  x(t) 
where  the  values  for  <| )i(L)  are  obtained  from  (6-18)  with  l  =  L. 


The  resulting  time  history  of  the  output  is  shown  below  in  Fig.  6-4.  Note  that  the  motion 
is  primarily  a  damped  sinusoid,  with  an  apparently  clean  harmonic  content.  The  power  spectral 
density  (PSD)  function  of  the  tip  displacement  is  shown  in  Fig.  6-5.  A  sampling  frequency  of 
600  Hz  was  used  in  the  simulation  so  the  response  of  all  modes  could  be  observed.  Note  that 
the  modal  frequencies  now  match  those  of  a  hinged-free  beam  because  of  the  introduction  of  the 
extra  degree  of  freedom  at  the  "fixed"  end,  i.e.  the  angle  0.  Mode  1  is  by  far  dominant,  being 
at  least  40  db  higher  than  mode  2.  Higher  modes  are  even  less  significant.  This  observation  will 
be  exploited  in  the  design  of  the  Kalman  filter  -  only  the  first  two  modes  will  be  modelled. 
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Simulated  tip  dlaplacarant^  2  cm  1C  Mode  1 


Fig.  6-4:  Simulated  tip  motion 


Fig.  6-5:  PSD  of  simulated  tip  motion 


6.3.3  Experimental  Corroboration 


The  model  of  the  tank  gun  barrel  was  partly  validated  by  a  series  of  field  trials  conducted 
on  a  Leopard  Cl  tank  with  a  105mm  gun  as  described  in  Bird  (1990).  Briefly,  the  barrel  was 
instrumented  with  a  number  of  accelerometers  and  gyroscopes  and  the  tank  was  driven  over 
various  types  of  natural  and  artificial  terrain.  Gyros  mounted  in  the  elevation  plane  at  the 
mantlet  and  the  muzzle  of  the  gun  provided  angular  rate  signals  that  were  sampled  at  60  Hz  and 
numerically  integrated  to  give  the  mantlet  and  muzzle  elevation  angles  such  as  those  shown  in 
Fig.  6-6  below.  An  estimate  of  the  elevation  flex  of  the  barrel  was  obtained  by  subtracting  the 
mantlet  angle  from  the  muzzle  angle  resulting  in  the  "derived  flex"  shown  in  Fig.  6-7.  The 
power  spectral  density  of  this  is  shown  in  Fig.  6-8.  Unfortunately,  the  signals  were  sampled  at 
only  60  Hz  so  PSD  information  is  only  valid  up  to  30  Hz.  The  dominant  mode  at  about  20  Hz 
appears  as  expected  but  higher  modes  are  cut  off  by  anti-aliasing  filters  at  30  Hz.  As  well,  the 
broadband  power  seen  in  the  1  Hz  to  5  Hz  region  is  that  of  the  tank  moving  over  the  ground  at 
moderate  speeds  and  was  not  included  in  the  simulations  of  the  last  section.  Nonetheless,  it  is 
encouraging  to  see  at  least  general  agreement  between  simulated  and  experimental  results.  In 
future  experiments  it  would  be  valuable  to  conduct  more  controlled  tests  with  the  tank  stationary 
with  barrel  motion  induced  by  initial  displacements  and  with  high  speed  data  recording  to  better 
validate  the  model. 

In  any  event,  it  will  be  assumed  that  the  model  is  a  reasonable  approximation  of  the 
dynamics.  Of  course,  a  number  of  assumptions  have  already  been  made  to  arrive  at  the  linear 
state  equations,  not  the  least  of  which  is  that  the  barrel  has  a  uniform  cross-section  along  its 
length.  If  this  assumption  is  not  made,  there  is  no  closed-form  solution  available  and  the 
dynamics  must  be  simulated  by  complex  numerical  techniques  such  as  finite  element  methods 
and  a  Kalman  filter  would  not  be  a  candidate  for  muzzle  pointing  prediction  algorithms.  The 
state  space  model  lends  itself  immediately  to  such  techniques. 
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Fig.  6-6:  Measured  barrel  angles 


Fig.  6-7:  Derived  barrel  flex 


Fig.  6-8:  PSD  of  flex 


6.4  A  KALMAN  FILTER  FOR  THE  DMRS 
6.4.1  Filter  Design 


This  section  describes  the  design  of  a  Kalman  filter  based  on  the  previously  developed  state 
equations.  The  filter  is  tested  on  its  ability  to  predict  the  muzzle  pointing  angle  of  the  real  data 
in  the  next  section. 

From  the  observations  of  the  last  section,  only  the  first  two  modes  of  vibration  are 
considered,  so  our  state  equation  is 


x  =  Fx  +  Gf 


where 


V 

*1 


/:= 


x 

0 

0 


and  where  the  coefficient  matrices  (for  the  2  mode  case)  work  out  to  be,  from  (6-21), 


F  matrix: 


-363.15  5 . 16462e-02 

1032.9  -0.14759 

164.9  -2.34511e-02 
1  0 
0  1 
0  0 


8.24476e-03 
-2 . 34511e-02 
-4 . 43573e-03 
0 
0 
1 


0 

0 

0 

0 

0 

0 


63551 
-1 . 81614e+05 
-28857 


0 

0 

0 


3 . 98487e+05 
-1 . 13344e+06 
-2 ,14389e+05 
0 
0 
0 


G  matrix: 


0.181574 

-0.516462 

-8.244756e-02 

0 

0 

0 


-0.516462 

1.47592 

0.234511 

0 

0 

0 


-8.244756e-02 
0.234511 
4 . 435732e-02 
0 
0 
0 


Define  the  product  Gf  as 


w  :=  Gf  = 


0.1816  X 
-0.5165  x 
-0.08244  x 
0 
0 
0 


(6-22) 


The  forcing  function,  x,  is  modelled  as  a  Gaussian  white  noise  stochastic  process  with  zero  mean 
(E  [t]  =  0)  and  power  spectral  density  (see  Appendix  B  for  a  justification  of  this): 
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E[T(/a)x(*6)J  =  5.0  Hta~tb) 


Thus  the  power  spectral  density  of  the  input  forcing  vector,  w,  is 


E[w(ta)wT(tb)] 


0.1816 

0.1816 

-0.5165 

-0.5165 

-0.08245 

0 

E[x(tayz(tb)] 

-0.08245 

0 

0 

0 

0 

0 

0.16  0  0  0  0  0 

0  1.3  0  0  0  0 

0  0  0.03  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 


8 (ta-tb)  :=  Q  Z(ta-tb) 


So  in  summary,  the  state  model  so  far  is 


x(t)  =  F  x(t)  +  w(t) 
E[w(i>T(riJ]  -  Q 


(6-23) 


This  is  exactly  the  form  of  the  continuous  time  stochastic  system  of  (2-2)  which  can  readily  be 
converted  to  the  discrete  time  system 


**♦1  =  ®xk  +  *k 


E[wkwJ]  =  Qk8kj 


according  to  the  discretization  procedure  (2-3). 


The  measurement  for  the  filter  will  be  the  angle,  y,  of  the  muzzle  with  respect  to  the 
mandet  (see  Fig.  6-9)  which  can  be  related  to  the  vertical  displacement  of  the  muzzle  from  the 
centreline  of  the  barrel,  T\(L,t),  and  the  length  of  the  gun,  L : 
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\|f(r)  =  sin-1 

L 


m  r\(L,t) 
L 


Fig.  6-9:  Measured  angle  from  muzzle  reference  system 

The  mantlet-mounted  laser  transceiver  of  the  muzzle  reference  system  measures  the  angular  flex 
of  the  barrel  using  a  reflecting  mirror  mounted  on  the  muzzle.  This  measurement  consists  of  the 
angle,  \jr,  plus  a  random  noise  component,  v*: 

z*  :«  ¥(**)  +  v* 

*  -  +  Vt 

L  * 

,  2 

-  '  Y.  oflXM  *  vt 

L  i«l 

=  i  [o  0  0  0  2  -2]  Xk  +  Vk 
L 

=:  H  xk  +  vk 

Finally  we  take  the  measurement  noise  to  be  zero  mean  and  covariance 
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E[v*Vy]  =  (0.2  mrad)2  bkj 
=  (0.0002  rad)2  5kj 
:=  Rk  by 

This  noise  level  was  chosen  based  on  expected  MRS  performance  and  also  happens  to  describe 
the  extent  of  the  noise  seen  in  the  field  data  of  the  next  section  quite  well 

6.4.2  Filter  Execution  on  Real  Data 


Now  that  the  specification  of  the  filter  is  complete,  the  derived  barrel  flex  from  the 
measured  field  trial  data.  Fig.  6-7,  was  used  as  the  data  sequence,  zk,  to  drive  the  filter  of  the 
previous  section.  It  is  difficult  to  evaluate  the  performance  of  the  filter  because  there  is  no  truth 
data  available  that  would  definitively  indicate  how  well  the  filter  is  predicting  the  future  muzzle 
pointing  angle.  All  we  have  is  the  sequence  of  noisy  measurement  data.  However  it  is  possible 
to  qualitatively  evaluate  the  filter  by  looking  at  the  predicted  muzzle  angle  relative  to  the 
measured  angle  as  a  function  of  time.  The  filter  was  used  to  predict  what  the  measured  flex 

would  be  3  steps  (or  50  ms,  since  the  measurement  rate  is  60  Hz)  ahead,  4*31*  =  Hk  4i*- 
(The  notation  here  is  the  same  as  in  Chapter  2.) 

Shown  in  Fig.  6-10  below  is  a  short  section  of  the  measured  barrel  flex  as  derived  from 
the  difference  of  the  two  gyros  at  opposite  ends  of  the  gun.  The  solid  line  is  the  actual 
measurement  sequence,  zk,  fed  to  the  filter  and  the  dotted  line  is  the  what  the  filter  predicted  the 
measurement  would  be  at  that  point  based  only  on  data  received  up  to  3  steps  (50  ms)  before  that 

3 

point,  i.e.,  fjyk_3  =  Hk_ 3  Ot_3  It  is  fairly  evident  from  Fig.  6-10  that  the  filter  is 

predicting  the  gun  flex  quite  well.  Given  that  the  filter  was  told  that  the  standard  deviation  of 
the  measurement  noise  was  0.2  mrad,  the  differences  between  measured  and  predicted  muzzle 
angles  all  appear  to  be  of  that  order,  as  could  be  expected,  in  that  section  of  the  data  sequence. 

However,  there  are  other  portions  of  the  measured  barrel  flex  data  sequence  that  were  not 
as  satisfactorily  predicted  by  the  filter.  One  of  these  is  shown  in  Fig.  6-11.  What  is  different 
about  this  portion  is  that  the  20  Hz  vibratory  mode  is  not  dominant  from  about  r=12.2  to  r=13.2. 
During  that  time  that  tank  likely  experienced  a  fairly  long  (1/2)  second  upward  acceleration 
followed  by  a  similar  downward  acceleration  that  essentially  biased  the  shape  of  the  gun  barrel 
slightly  while  the  accelerations  were  occurring.  This  effect  was  not  included  in  the  barrel  model 
used  by  the  Kalman  filter  and  hence  the  predicted  measurements  during  these  times  are  not 
particularly  good.  In  fact  they  tend  to  lag  behind  by  3  or  4  samples.  After  the  acceleration 
transients  die  away  at  about  13.2  sec,  the  filter  returns  to  fairly  good  prediction  performance. 
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Measured  and  3  step  ahead  (SO  ms)  filter-predicted  flex 


t(s) 

Fig.  6*10:  Real  baiTel  flex  measurements  and  3-step  ahead  Kalman  filter 

prediction.  This  shows  fairly  good  muzzle  angle  prediction  when 
the  20  Hz  mode  is  the  dominant. 


Measured  and  3  step  ahead  (50  ms)  filter-predicted  flex 


t(s) 

Fig.  6-11:  Real  barrel  flex  measurements  and  3-step  ahead  Kalman  filter 

prediction.  This  shows  poorer  muzzle  angle  prediction  when 
unmodelled  accelerations  dominate. 
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6.4.3  Discussion  of  Filter  Results 

The  Kalman  filter  design,  as  presented,  is  not  yet  of  sufficient  fidelity  to  include  in  a  real¬ 
time  predictive  DMRS  since  it  does  not  accurately  model  the  low  frequency  effects  caused  by 
sustained  accelerations  as  manifested  in  Fig.  6-11.  This  is  not  an  insurmountable  problem. 
Relatively  modest  changes  to  the  state  equations,  process  noise  models,  and  software  simulator 
and  careful  tuning  of  the  noise  parameters  would  yield  a  filter  that  could  predict  the  muzzle  angle 
through  these  sustained  accelerations  as  accurately  as  the  current  model  can  predict  the  higher 
frequency  vibrations.  This  has  not  yet  been  done  since  the  filter  so  far  has  been  designed  around 
the  angular  flex  data  derived  from  gyroscopes,  while  the  DMRS  system  that  will  be  installed  on 
the  tank  will  be  laser  and  mirror  based  This  is  such  a  sufficient  difference  that  fine-tuning  the 
model  on  the  existing  data  may  not  yield  a  suitable  filter  for  data  from  the  laser  system.  In  any 
event,  the  results  so  far  have  been  quite  encouraging  and  the  modelling  efforts  in  this  chapter 
should  be  a  good  launching  point  for  further  studies  on  predictive  muzzle  referencing  systems. 
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7.0  A  POTENTIAL  ADM  ARCHITECTURE 


This  chapter  presents  a  preliminary  proposal  for  an  overall  architecture  of  an  Advanced 
Development  Model  for  a  highly  advanced,  fire-on-the-move  fire  control  system  for  future  direct- 
fire  support  vehicles.  The  previous  chapters  have  dealt  in  detail  with  some  of  the  more  pressing 
problems  such  as  trajectory  prediction  and  dynamic  muzzle  referencing.  However  it  may  be 
useful  to  briefly  step  back  and  look  at  the  overall  system  as  a  whole. 

Fig.  7-1  is  an  attempt  at  showing  both  the  major  system  components  and  the  sequence  of 
events  that  would  occur  as  the  data  travels  through  the  system  from  the  time  of  initial  target 
detection  until  the  weapon  is  fired. 

First,  we  assume  that  there  will  be  two  operators  of  the  system.  In  current  generations  of 
vehicles,  they  are  called  Commander  and  Gunner  but  we  shall  use  the  more  generic  terms 
Observer  and  Engager.  The  diagram  has  a  top  portion  for  the  equipment  associated  with  the 
observer’s  duties  and  a  bottom  portion  for  equipment  associated  with  the  engager’s.  Each 
operator  has  his  own  imaging  sensor(s)  that  can  be  aimed  independently  of  the  weapon  and  each 
other.  Associated  with  each  operator  station,  is  a  high-speed,  dedicated  image  processing 
computer  that  can  be  configured  by  the  operator  to  aid  in  image  enhancement  (by  contrast 
stretching,  motion  highlighting,  etc.)  image  stabilization  (jitter  removal  by  registering  the  image 
against  a  stationary  background),  target  extraction  (by  highlighting  potential  targets  and  extracting 
the  more  promising  ones)  and  target  tracking  (by  running  Kalman  filters  similar  to  those  in 
Chapter  4)  to  assist  in  finding  the  target  in  the  next  scan. 

A  sensor  control  computer  will  be  used  to  control  the  functions  of  each  image  processor 
and  act  as  the  human  interface  to  the  system.  This  control  computer  can  be  used  to  suggest  a 
prioritization  of  the  targets  to  the  observer,  possibly  assist  in  their  identification,  accept  all 
operator  hand  or  voice  inputs,  and  cue  the  imaging  sensor  to  follow  a  specific  target  or  the 
operator’s  hand  control.  The  control  computers  of  the  two  stations  would  be  in  constant 
communication  so  that  targets  could  be  handed  off  from  the  observer  to  the  engager,  or  in 
emergencies  could  redirect  the  output  of  one  image  processor  to  the  other  operator’s  display. 

When  it  is  decided  that  a  particular  target  is  to  be  engaged,  the  engager’s  control  computer 
would  concentrate  solely  on  that  target.  It  would  run  a  more  sophisticated  Kalman  tracking  filter, 
perhaps  with  more  states,  maneuver  detectors,  etc,  to  attempt  to  estimate  very  accurately  the  true 
state  of  the  target  It  would  drive  the  imaging  sensor  to  closely  follow  that  target.  An  interesting 
possibility  is  that  it  might  be  able  to  estimate  the  probability  of  hitting  the  target  based  on  its 
estimate  of  the  target  acceleration,  for  example. 

Provided  the  engager  is  satisfied  with  his  chance  of  success,  the  laser  range  finder  (LRF), 
environment  sensors  (wind,  air  pressure,  temperature,  etc)  would  be  activated  by  the  control 
computer  and  the  ballistic  solution  would  be  calculated  based  on  the  data  from  these  sensors  and 
the  target  velocity/acceleration  available  from  the  target  state  estimator.  The  gun  control  servos 
would  slew  the  gun  to  the  proper  lead  angles  in  preparation  for  firing.  At  this  point,  an  optional 
interrogation  of  the  target  might  take  place  based  on  an  optical  IFF  (Identify  Friend  or  Foe) 
system  and  the  computer  would  inform  the  engager  that  firing  can  commence. 
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Fig.  7-1:  A  Potential  ADM  Architecture 


When  the  engager  presses  the  firing  button,  there  may  be  a  slight  delay  as  the  barrel  flex 
data  from  the  DMRS  (Dynamic  Muzzle  Reference  System)  laser  on  the  gun  barrel  (and  possibly 
other  vehicle  dynamics  from  the  onboard  Inertial  Navigation  System  (INS))  is  processed  through 
a  predictive  Kalman  filter  similar  to  that  of  Chapter  6  and  a  barrel  zero  crossing  time  is 
predicted.  The  charge  would  be  ignited  at  a  time  such  that  the  shell  could  be  expected  to  exit 
the  muzzle  at  a  zero  crossing  time. 

After  the  results  of  the  initial  shot  are  observed,  the  observer  will  have  the  next  target  ready 
for  the  engager.  The  target  is  handed  off  to  the  engager’s  control  computer  and  the  process 
repeats. 

From  this  very  cursory  design,  it  can  be  seen  that  such  a  system  is  very  flexible,  powerful 
and  very  comprehensive.  The  algorithms  discussed  in  this  report  make  up  just  two  of  the  overall 
components,  the  optimal  firing  time  solution  and  the  target  state  estimation.  Much  more  work 
remains  to  be  done  oefore  such  a  system  is  ready  to  be  fielded. 


74 


8.0  SUMMARY 


This  report  has  summarized  some  applications  of  Kalman  filters  in  the  fire  control  systems 
of  future  direct-fire  support  vehicles,  especially  in  the  high  dynamic  moving-target/moving- 
platform  scenario.  Two  areas  were  studied  in  some  detail,  target  state  estimation  and  dynamic 
muzzle  referencing.  A  number  of  modelling  studies  were  conducted,  simulation  software  was 
written  and  exercised,  candidate  filters  were  designed  and  implemented  and,  where  possible,  real 
data  was  used  to  test  the  resulting  filters.  The  results  show  sufficient  potential  to  conclude  that 
such  techniques  can  be  extremely  useful  in  predicting  the  dynamics  of  both  the  target  and  the 
firing  platform.  An  overall  architecture  of  an  advanced  development  model  was  proposed  to 
establish  the  context  in  which  the  filters  of  this  report  can  play  a  vital  role. 
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APPENDIX  A.  DERIVATION  OF  Q  VALUES  FOR  THE  TRACKING  FILTERS 


The  general  continuous-tune  model  of  target  motions,  from  (2-2),  is 

xc(t)  =  Fxc(t)  +  wc(t) 


(A-l) 


with  the  covariance  of  the  driving  process  given  by 

E  {wc(t)wc(x)T)  =  Qc(f)6(f-t) 

If  we  look  only  at  the  velocity  state  in  one  horizontal  direction,  say  vYG,  the  randon  walk  velocity 
(RWV-r)  continuous  time  velocity  model  is  driven  by  white  noise  (from  Eqn.  (3-14)): 

Vyc(0  =  waYG(t)  (A-2) 


We  wish  to  estimate  a  value,  Qc,  for  the  covariance  of  w  .  We  will  approach  the  problem  first 

from  the  discrete-time  point  of  view  and  then  return  to  continuous  time.  The  discrete  time 
version  of  (A-l)  is,  as  shown  by  (2-1): 

**♦1  «  xk  +  wk  <A-3> 

where,  from  (2-3), 

<pk  =  (A-4) 

But  since  F= 0  for  the  velocity  state  (A-2)  of  the  RWV  model,  then  d>  =  1.  The  discrete-time 
velocity  model  is  thus 


vYc(lk*0  =  VYC^0  +  wk  (A*5) 

where  we  have  simplified  notation  by  redefining  wk  as  the  scalar  noise  process  driving  the  vYG 
velocity  state.  The  covariance  of  this  discrete  driving  noise  is  also  given  by  (2-3): 

E  {w*}  =  Qk  =  j[^+1ef(‘t*I'X)Gc(x)eFr(‘**,‘/t)dx  (A-6) 

but  since  F= 0,  this  simplifies  to 
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Qk  =  (! k*lQc«W 

Jtk 


(A-7) 


which,  for  sufficiently  small  time  steps  (tk+1  -  tk  -»  0),  can  be  approximated  by 

Qk  58  2C(**)  (**+i  -  **) 


(A-8) 


r 


Thus  by  estimating  a  value  for  Qk  =  Efw*2},  we  can  obtain  the  Qc  we  are  looking  for. 

Let  us  assume  that  a  reasonably  maneuverable  ground  target  is  to  be  tracked,  say  one  that 
can  accelerate  from  0  to  100  km/hr  in  5  seconds,  or  equivalently,  can  execute  a  90  degree  turn 
at  100  km/hr  in  5  seconds  (both  of  which  are  quite  realistic  for  modem  fighting  vehicles).  This 
is  equivalent  to  an  acceleration  of  5.55  m/s2.  With  a  discretization  time  step  of  tk+1  -  tk  -  0.1s, 
this  corresponds  to  a  velocity  change  of  0.555  m/s  in  0.1s.  Thus  we  will  model  the  noise  driving 
the  discrete  velocity  state  as  an  uncorrelated  Gaussian  sequence  with  zero  mean  and  covariance 


Qk  =  (0.555  m/s)2 
=  0.3  m2/s2 


(A-9) 


so  that  from  (A-8), 


Qc  ~  Qk  !  (r*+i 

=  0.3  /  0.1  (A- 10) 

=  3  m2/s3 


In  the  vertical  direction,  we  assume  the  vehicle  is  only  a  tenth  as  maneuverable  as  it  is 
horizontally.  Thus  Qk  in  the  vertical  is  (0.0555  m/s)2  so  the  corresponding  Qc  is  0.03  m2/s3. 
These  are  the  values  used  in  the  Qc  matrix  in  the  RWV(r)  simulations  of  Chapters  4  and  5. 

For  the  angular  rate  model,  RWV(s),  a  similar  derivation  is  used.  Since  the  range  to  the 
target  is  unknown,  however,  a  nominal  range  of  about  2000m  is  used  to  derive  the  Q  values. 
(This  is  not  atypical  for  land  engagements.)  A  cross-range  linear  velocity  change  of  0.555  m/s 
over  0.1s  at  a  range  of  2000m  results  in  an  angular  velocity  change  of  about  0.25  mrad/s.  Thus 
the  discrete  noise  driving  the  horizontal  angular  velocity  state  will  be  assumed  to  have  a 
covariance  of 


Qk  =  (0.00025  rad/s)2 
=  6.25  x  10"8  rad2/s2 


(A-ll) 
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so  that  again  from  (A-8), 


Qc  "  Qk  !  (*k+ 1  “**) 

=  6.25  x  10-8  /  0.1  <A*12> 

=  6.25  x  10“7  rad2/s3 

The  corresponding  vertical  angular  velocity  covariance,  assuming  the  similar  1/10  maneuver¬ 
ability,  results  in  a  Qc  of  6.25  x  10*9  rad2/s3. 


■» 
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APPENDIX  B.  DERIVATION  OF  Q  VALUES  FOR  THE  MRS  FILTER 


To  estimate  the  Q  matrices  for  the  muzzle  reference  filter,  we  take  a  similar  approach  as 
for  the  tracking  filters  in  Appendix  A.  First  we  must  determine  a  level  a  "maneuverability"  of 
the  gun  barrel  by  looking  at  some  of  the  collected  field  data.  For  the  same  sequence  of  data 
shown  in  Fig.  6-6,  the  angular  rate  and  acceleration  of  the  mantlet  are  shown  in  Fig.  B-l  and 
Fig.  B-2  respectively. 


» 


I 


< 


Fig.  B-l:  Mantlet  angular  rate 


Fig.  B-2:  Mantlet  angular  acceleration 


From  Fig.  B-2,  a  not  unreasonable  assumption  for  the  level  of  maneuverability  of  the  mantlet  can 
be  rather  arbitrarily  chosen  as  3  rad/s2.  For  the  60  Hz  sampling  rate  of  the  data  in  question,  this 
corresponds  to  an  angular  velocity  change  of  about  0.051  rad/s  in  1/60*  of  a  second.  Thus  we 
will  model  the  noise  driving  the  discrete  angular  velocity  state  as  an  uncorrelated  Gaussian 
sequence  with  zero  mean  and  covariance 


Qk  =  (0.051  rad/s)2 
=  0.0026  rad2/s2 


(B-l) 


so  that  from  (A-8), 


Qc  **  2*  /  (**+i 

=  0.0026  /  0.0166667  (B-2) 

=  0.16  rad2/s3 

From  the  first  element  of  the  vector  equation  relating  the  noise  driving  the  angular  velocity  state 
with  the  noise  on  the  applied  torque  x,  Eqn.  (6-22),  we  can  conclude 
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Ema)x(tb))  =  0.16  /  0.18162  6 (ta  -  tb) 
=  5  5  (ta  -  tb) 


(B-3) 


This  value  is  then  used  to  derive  the  remainder  of  the  Q  elements  of  Chapter  6. 
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